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PREFACE 

The  work  reported  here  was  performed  by  the  Optical  Science 
Laboratory  of  the  Advanced  Concepts  Division,  Environmental  Research 
Institute  of  Michigan  (ERIM).  The  work  was  sponsored  by  the  Office  of 
Naval  Research  and  the  Naval  Research  Laboratory  under  Contract 
N00014-86-C-0587 .  At  NRL,  the  technical  monitor  was  Dr.  Jerry  A. 
Blodgett. 

This  final  technical  report  covers  work  performed  from  1  August 
1987  to  31  January  1989.  The  principal  investigator  at  ERIM  was  Jack 
Cederquist.  Major  contributors  to  this  work  were  Jack  Cederquist, 
James  R.  Fienup,  Ann  M.  Kowalczyk,  Joseph  C.  Marrnn.  and  Kirk  S. 
Schroeder. 
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INTRODUCTION 


1.1  BACKGROUND 

The  ability  to  discriminate  between  targets  (e.g.,  re-entry 
vehicles)  and  decoys  is  of  great  importance  to  our  nation's  strategic 
defense.  One  approach  for  accomplishing  such  discrimination  is  to 
obtain  a  fine-resolution  image.  Conventional  approaches  for  fine- 
resolution  optical  (UV  to  near  IR)  imaging  of  objects  at  great 
distances  require  large-diameter  receiver  optics  with  near-diffraction- 
limited  imaging  performance.  A  very  stable  receiver  support  structure 
and/or  adaptive  correction  of  optical  misalignment  to  within  a  fraction 
of  a  wavelength  is  also  required  to  maintain  this  near-diffraction- 
limited  performance.  Further,  if  the  receiver  must  operate  over  a 
large  f ield-of-regard,  prohibitively  complex  hardware  is  required.  In 
our  research  program  we  have  investigated  a  meaningful  and  realizable 
alternative  solution  to  meet  the  challenge  of  fine-resolution  imaging. 

For  coherently  illuminated  objects,  a  class  of  alternative, 
unconventional  approaches  to  fine-resol ution  imaging  uses  a  large  array 
of  small  diameter  receivers  to  make  Fourier  intensity  measurements  in 
the  receiver  aperture.  From  these  measurements  and  additional  low- 
resolution  imagery  and/or  a  priori  information  about  the  object,  it  is 
possible  to  compute  a  fine-resolution  image  using  a  phase  retrieval 
algorithm.  The  intensity  measurements  can  be  made  for  objects  over  a 
large  fiel d-of-regard.  In  comparison  to  conventional  approaches,  this 
method  greatly  reduces  the  receiver  hardware  requirements  in  exchange 
for  increasing  the  data  processing  necessary  to  compute  the  image.  In 
addition,  certain  object  parameters,  such  as  rotation  rate,  can  be 
determined  from  the  Fourier  intensity  measurements  without  forming  an 
image. 
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1.2  OVERVIEW  OF  ACCOMPLISHMENTS 

This  report  describes  the  results  of  an  18  month  program  for 
development  of  unconventional  approaches  to  image  reconstruction  and 
parameter  estimation.  In  this  section  the  principal  results  of  the 
research  program  will  be  briefly  summarized.  They  are  reported  in 
detail  in  the  sections  that  follow.  Sections  2  and  3  give  research 
results  in  the  area  of  parameter  estimation  and  Sections  4  and  5 
describe  image  reconstruction  research.  Section  6  reports  on  the 
investigation  cf  system  requirements. 

The  research  conducted  in  the  area  of  parameter  estimation  concerns 
two  types  of  object  parameters.  The  first,  which  is  described  in 
Section  2  is  object  rotation  rate  and  the  second,  object  separation,  is 
the  subject  of  Section  3.  In  the  area  of  rotation  measurement  we  have 
developed  a  simple  model  that  allows  us  to  compute  the  space-time 
correlation  function  of  the  speckle  patterns  from  objects  of  arbitrary 
shape  and  surface  material.  This  computational  model  was  confirmed 
experimentally  on  cylinders  with  a  variety  of  surface  materials.  We 
also  demonstrated,  both  in  theory  and  through  experiments,  rotation 
rate  measurements  of  multiple  objects  illuminated  simultaneously. 
Other  topics  covered  include  the  robustness  of  rotation  measurement  at 
low  light  levels  and  a  comparison  of  the  speckle  methods  for  rotation 
measurement  with  heterodyne  methods.  Measurement  of  object  separation 
rate  from  the  Fourier  transform  of  speckle  pattern  data  is  described  in 
Section  3.  Both  a  theoretical  treatment  of  separation  measurement  from 
speckle  intensity  and  experimental  results  that  demonstrate  the  method 
are  given. 

Ordinarily  it  is  very  difficult  to  reconstruct  an  image  of  a 
complex-valued  object  from  the  modulus  of  its  Fourier  transform  ( i . e . , 
retrieve  the  Fourier  phase)  except  for  some  special  cases.  In  Section 
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4,  a  two-step  approach  is  described  for  reconstructing  high-resolution 
images  of  a  general  object  from  Fourier  intensity  data  using 
additionally  a  low-resolution  intensity  image.  First,  the  Fourier 
phase  over  the  small  aperture  is  retrieved  using  the  Gerchberg-Saxton 
algorithm.  Then  that  phase  is  used,  in  conjunction  with  the  Fourier 
modulus  data  over  a  large  aperture  together  with  a  support  constraint 
on  the  object,  to  reconstruct  a  fine-resolution  image  (retrieve  the 
phase  over  the  large  aperture)  by  the  iterative  Fourier  transform 
algorithm.  A  series  of  simulations  that  demonstrate  image 
reconstruction  and  test  the  sensitivity  of  the  Gerchberg-Saxton 
algorithm  to  photon  noise  are  described.  Experiments  were  also 
conducted  to  demonstrate  image  reconstruction  from  Fourier  intensity 
data  obtained  in  laboratory  experiments.  The  details  of  the 
experimental  setup,  system  calibration,  and  successful  reconstruction 
results  for  a  simple  two  part  object  are  given  in  Section  5.  Finally, 
Section  6  of  this  report  summarizes  work  done  to  specify  the  laser, 
detector,  and  data  processor  requirements  of  a  deployed  system  for 
image  reconstruction  and  parameter  estimation. 

1.3  SUMMARY  REMARKS  AND  RECOMMENDATION 

In  summary,  we  have  succeeded  in  developing  important  new  image 
reconstruction  techniques  and  speckle-based  target  parameter 
measurement  techniques.  These  techniques  have  a  strong  potential  for 
use  in  the  SDI  midcourse  discrimination  application.  Further,  it  is 
important  to  note  that  these  techniques  do  not  require  the  use  of  large 
diffraction-limited  optics  or  heterodyne  detection  and  should  therefore 
lead  to  more  reliable,  lower  cost,  and  lighter  weight  system  designs. 
We  recommend  additional  research  to  further  develop  and  apply  these 
techniques . 
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MEASUREMENT  OF  OBJECT  ROTATION  USING  LASER  SPECKLE 

2.1  INTRODUCTION 

Consider  the  speckle  pattern  produced  by  a  fully  illuminated,  3-D 
object  as  shown  in  Fig.  2-1.  It  is  assumed  that  the  surface  of  the 
object  is  sufficiently  rough  to  produce  a  fully  developed  speckle 
pattern  in  the  observation  plane.  When  the  object  rotates  about  a  fixed 
axis,  the  dominant  effect  observed  is  that  the  speckle  pattern 
translates  perpendicular  to  the  object's  axis  of  rotation  and  that  the 
distance  of  speckle  translation  is  proportional  to  the  amount  of  object 
rotation.  Finer  observation  reveals  that,  along  with  the  translation, 
the  speckle  pattern  decorrelates  or  'boils.'  The  degree  of  boiling  is 
primarily  dependent  on  the  amount  of  object  rotation  and  factors  such  as 
the  locations  of  the  source  and  observation  plane  as  well  as  the 
underlying  shape  of  the  object. 

In  this  section  we  analyze  a  simple  computational  model  useful  for 
calculating  space-time  correlation  functions  for  speckle  from  rough 
rotation  objects.  These  theoretically  calculated  correlation  functions 
can  be  used  to  assess  the  ability  to  make  remote  measurements  of  an 
object's  rotation  rate  and  shape  from  speckle  intensity  correlation. 

Speckle  from  rotating  objects  in  free-propagation  arrangements  has 
been  investigated  by  several  authors.  George  [2.1]  has  provided  a 
detailed  analysis  of  the  space-time  correlation  properties  of  speckle  in 
the  far  field  of  the  object.  He  considered  point-source  illumination  of 
the  entire  object  and  calculated  the  space-time  correlation  function. 
An  analysis  of  the  Doppler  spectrum  obtained  by  interfering  the 
scattered  light  with  a  local  oscillator  was  also  given. 
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Figure  2-1.  Example  of  an  optical  system  used  to  observe  dynami 
speckle  from  rough,  rotating  objects. 


Erdmann  and  Gellert  [2.2]  and  later  Takai  et  al.  [2.3]  analyzed 
speckle  from  rotating  objects  illuminated  by  a  laser  spot.  Erdmann  and 
Gellert  developed  an  expression  for  the  space-time  intensity  correlation 
function  in  the  far  field.  Takai  and  co-workers  considered  a  single 
detector  In  the  near  field  and  derived  an  expression  for  the  temporal 
intensity  correlation  function.  Both  groups  presented  experimental 
measurements  to  verify  the  theoretical  results. 

Leader  [2.4]  has  investigated  temporal  properties  of  the  speckle  by 
analyzing  the  frequency  spectrum  of  the  detected  light.  Smith  [2.5] 
developed  criteria  to  indicate  when  Doppler  frequency  shifts  of  the 
scattered  light  dominate  the  speckle  effects. 

Hayashi  and  Kitagawa  [2.5]  reported  a  technique  to  measure  the 
rotation  angle  of  a  cylinder  using  near-field  speckle.  They  derived  a 
relationship  in  which  the  speckle  translation  is  equal  to  the  object 
rotation  multiplied  by  a  constant  and  confirmed  their  results  with 
experiments. 

From  the  investigations  mentioned  above  it  becomes  apparent  that 
calculations  of  the  space-time  correlation  function  of  the  optical  field 
can  be  difficult,  particularly  when  the  object  has  a  complicated 
underlying  shape.  In  this  section  we  present  a  relatively  simple  ray- 
trace  method  for  computing  space-time  correlation  functions  of  dynamic 
speckle.  Here  we  consider  only  rotation  about  a  fixed  axis.  However, 
the  method  can  be  applied  to  more  complicated  motions.  The  theoretical 
basis  for  the  ray-trace  method  is  presented  in  Section  2.2.  In  Section 
2.3  the  method  is  used  to  compute  the  space-time  correlation  function 
for  a  series  of  2-D  disk  objects;  these  calculations  allow  one  to 
develop  an  Intuitive  understanding  of  the  amounts  of  speckle  boiling  and 
translation  to  expect  for  general  objects.  In  Section  2.4  experimental 
calculations  of  the  correlation  function  for  speckle  from  a  rotating 
cylinder  are  presented  and  compared  to  theory.  In  Section  2.5  we 
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consider  multiple  objects  rotating  at  different  rates  and  Section  2.6 
contains  comments  on  the  application  of  this  model  to  objects  with 
complicated  underlying  shape.  Section  2.7  contains  analysis  of  the  SNR 
for  rotation  measurement  at  low  light  level  and  in  Section  2.8  the 
intensity  detection  methods  presented  in  this  paper  are  compared  to 
heterodyne  methods  particularly  at  low  light  level. 

2.2  ANALYSIS  OF  SPACE-TIME  CORRELATION 

A  single  polarization  component  of  the  optical  field  scattered  by  a 
rough  object  can  always  be  written  as  a  sum  of  contributions  from 
discrete  scattering  cells  on  the  surface.  The  optical  field  at  an 
observation  point  p  is  thus 

U(p)  =  X]  exp[i>k]  (2-1) 

where  lA^I  and  are  the  amplitude  and  phase  of  the  contribution  from 
the  kth  scattering  cell  respectively.  The  value  of  lA^I  is  determined 
by  the  scattering  cross-section  (or  inclination  factor)  of  the 
particular  cell  and  the  strength  and  uniformity  of  the  illuminating 
field.  The  phase,  ^ ,  is  given  by  the  optical  path  length  of  the  light 
as  it  travels  from  the  source  to  the  scattering  cell  to  the  detection 
point.  In  this  treatment  we  ignore  phase  shifts  caused  by  reflection 
from  the  surface  and  if  the  object  depolarizes  the  incident  light  it  is 
assumed  that  a  polarization  analyzer  is  placed  in  front  of  the  detector. 

To  analyze  the  dynamics  of  the  speckle  pattern  produced  by  target 
rotation  we  will  use  a  quasi-static  analysis  of  the  optical  field  and 
thus  concentrate  on  classical  effects  caused  by  position  changes  of  the 
scattering  cells.  It  is  important  to  note  that  speckle  dynamics 
predicted  by  this  quasi-static  analysis  can  equivalently  be  interpreted 
as  self-doppler  or  autodyne  effects  [2. 7,2. 8]. 
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From  Eq.  (2-1)  we  see  that  changes  in  the  optical  field  caused  by 
object  rotation  result  from  the  following  three  effects; 

a)  The  particular  scattering  cells  that  enter  the  summation  in  Eq. 
(2-1)  change  as  scatterers  move  into  and  out  of  the  illuminated 
region  of  the  object.  The  severity  of  this  effect  is  dependent  both 
on  the  extent  of  the  illuminated  region  and  shadowing  caused  by  the 
underlying  shape  of  the  object.  Shadowing  caused  by  the  surface 
roughness  must  also  be  considered  at  large  bistatic  angles. 

b)  The  amplitude  of  a  scatterer  lA^I  can  change  because  of  angular 
dependence  of  the  cross-section.  Consider,  for  example,  a  faceted 
surface.  The  cross-section  is  very  large  when  a  facet  normal 
directs  incident  rays  from  the  source  to  the  detection  point  and  can 
drop  off  greatly  as  the  facet  normal  rotates. 

c)  The  phases  of  the  contributions  change  because  of  motion  of  the 
scattering  cells  caused  by  object  rotation. 

Our  analysis,  in  both  theory  and  experiment,  has  shown  that,  for 
most  objects,  speckle  boiling  and  translation  can  be  accounted  for  by 
the  phase  changes  described  In  Item  c)  alone.  This  conclusion  is  also 
supported  by  Leader  [2.4],  The  other  effects  can  be  significant  if  the 
underlying  shape  of  the  object  Is  very  complicated,  the  illuminated 
region  of  the  object  is  small,  or  the  cross  section  of  the  scattering 
cells  has  significant  angular  dependence.  However,  for  the  remainder  of 
this  analysis  we  concentrate  on  speckle  dynamics  that  result  from  phase 
changes  and  proceed  with  the  assumption  that  the  other  effects  are 
negligible.  The  validity  of  this  assumption  is  tested  by  the 
experiments  presented  in  Section  2.4. 
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To  quantify  the  dynamics  of  the  speckle  pattern  the  space-time 
correlation  function  of  the  optical  field  is  used.  Following  the 
discussion  above,  our  expression  for  the  correlation  function  is 
dependent  on  object  rotation  only  through  the  phase  of  the  scattering 
cells.  The  space-time  correlation  function  of  the  optical  field  is  then 

<U,  U*>  =  H  H  <IAkl  lAml  exp[i(^kl  -  ^m2)]>  ,  (2-2) 

km 

where  the  subscripts  1  and  2  denote  space  time  coordinates  (p^,t^)  and 
(p2 r t~)  respectively.  The  phase  term  <p^  represents  the  phase  of  the 
light  from  propagation  from  the  source  to  the  k^  scattering  cell  to  the 
detection  point  p^  at  time  t1  before  object  rotation.  The  phase  is 
defined  analogously. 

To  proceed  with  this  calculation  we  follow  Goodman  [2.9]  and  place 
tht  following  requirements  on  the  field  contributions: 

1)  The  amplitude  and  phase  of  the  kth  scatterer  are  statistically 
independent  of  each  other  and  of  the  amplitudes  and  phases  of  all 
other  scatterers. 

ii)  The  phases  are  uniformly  distributed  over  the  fundamental  interval 

[-».*] • 


With  these  requirements  satisfied  we  adopt  the  delta  correlation  model, 

<'V  'V  “Pf’t'kl  -  ■  5kn,  <'V2>  ■  <2-3> 

where  <5km  is  the  Kronecker  delta  function.  The  space-time  correlation 
function  of  the  optical  field  then  reduces  to 
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where  we  have  substituted  A^k  =  ^  ^ . 

To  calculate  the  correlation  function  we  will  modify  the  sum  in  Eq. 
(2-4);  instead  of  summing  over  the  scattering  cells  indexed  by  k,  we 
will  sample  the  expected  cross-section,  <IA^I  >,  and  phase  change,  A f>t  , 
over  a  discrete  grid  of  points  and  sum  over  these  samples.  This  new  sum 
is  equivalent  to  tracing  a  series  of  rays  that  coarsely  sample  the 
object  and  adding  up  the  corresponding  magnitudes  and  phases;  for  this 
reason  we  refer  to  the  technique  as  a  ray-trace  method.  Using  modified 
notation,  the  correlation  function  as  calculated  using  the  ray-trace 
method  is 


<u!  u2>  =  H  <IA(xk)l2>  exp[iA^(xk)]  (2-5) 
k 

where  xk  denotes  the  ray  coordinates.  There  are  several  subtle  issues 
associated  with  choosing  the  proper  ray  coordinates  for  sampling  a  given 
object,  tor  this  investigation  we  have  used  a  sampling  grid  with 
equally  spaced  ray  coordinates;  the  corresponding  grid  for  a  2-D  object 
is  shown  in  Fig.  2-2  (the  extension  to  3-D  objects  is  straightforward 
and  will  not  be  discussed  here).  The  line  of  the  sampling  grid  passes 
through  the  axis  of  rotation  of  the  object  and  is  perpendicular  to  the 
vector  connecting  the  axis  of  rotation  and  detection  point  Pj.  The 
sampling  rays  are  traced  from  the  source  to  points  on  the  surface  with  x 
axis  positions  given  by 


xk  =  xQ  +  k  Ax 


(2-6) 
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Figure  2- 


.  Sampling  grid  of  equally  spaced  ray  coordinates  used  to 
compute  the  speckle  correlation  function.  The 
relationship  of  this  grid  to  the  object  and  detection 
points  is  shown. 
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where  xQ  designates  the  beginning  of  the  object  as  shown  in  Fig.  2-2  and 
Ax  is  the  sample  spacing.  The  axis  of  rotation  is  located  at  x  =  0.  A 
computer  program  used  to  calculate  the  correlation  function  would  first 
determine  the  coordinate  yk  on  the  surface  that  corresponds  to  the 
coordinate  xk  as  shown  in  Fig.  2-2.  The  program  would  then  determine 
the  cross-section,  < I A (x^) 1 2>  of  the  object  at  (xk,yk)  and  compute  the 
value  of  the  phase  corresponding  to  the  optical  path  of  a  ray  that 
travels  from  the  source  to  (xk,yk)  and  to  the  detection  point  p.. 
Rotation  of  the  object  is  then  accounted  for  by  rotating  the  coordinates 
°f  (xk,y.)  as  follows 


xk  =  xk  cos  A0  +  yk  sin  A 9  (2-7) 

yk  =  -xk  sin  A0  +  yk  cos  A 9  ,  (2-8) 

where  the  object  rotates  about  the  point  x  =  y  =  0  and  A0  is  the  amount 
of  object  rotation  given  by  A0  =  fl(t2  -  t^  with  fl  being  the  angular 
rotation  rate.  More  complicated  object  motions  can  be  handled  by 
modifying  the  transformation  in  Eqs.  (2-7)  and  (2-8).  The  value  of 
A^(xk)  is  found  by  computing  the  difference  (^kl  -  ^k2)  where  ^k2  is  the 
phase  of  the  optical  path  between  the  source,  the  scatterer  at  (xk,  yk) 
and  the  detection  point  p2. 

With  knowledge  of  < I A (xk) I  >  and  A^(xk)  at  many  points  on  the  object 
we  can  compute  the  correlation  function  of  the  optical  field  using  Eq. 
(2-5).  However,  it  is  often  more  convenient  to  use  the  normalized  form 
of  the  correlation  function  which  is  given  by  [2.9] 

<U,  u2> 

P 12  _  r  *  *11/2"  (2-9) 

12  [<U1  UT>  <U2  U2>]1/2 
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Values  of  fall  in  the  range  0  S  1^1  i  1  with  l/i^l  =  1 

corresponding  to  pure  speckle  translation  and  l/i^l  <  1  indicating 
boiling.  After  substituting  Eq.  (2-5),  is  given  by 


X!  <IA(xk)l2>  exp[iA^(xk)] 
k  _ _ 

H  O  A(x.  )  I  2> 
k  K 


(2-10) 


From  Eq.  (2-10)  it  is  apparent  that  pure  speckle  translation  results 

when  A^(xk)  is  a  constant  for  all  xk  and  that  the  amount  of  speckle 

decorrelation  is  dependent  on  the  deviation  of  A^(xk)  from  being  a 

constant.  There  is  strong  similarity  between  A^(x, )  and  wavefront 

2  K 

aberration  of  imaging  systems.  In  fact  l/j^l  is  analogous  to  the 
normalized  intensity  from  the  diffraction  theory  of  aberrations  [2.10]. 
The  normalized  intensity  is  used  to  quantify  the  aberration  content  of 
an  imaging  system. 

In  summary,  we  have  developed  a  model  in  which  speckle  dynamics  can 
be  accounted  for  by  phase  changes  of  the  scattered  light  caused  by 
object  rotation.  Speckle  translation  results  when  the  phase  changes 
imparted  by  object  rotation  are  cancelled  by  the  phase  change  caused  by 
translation  of  the  detection  point.  Speckle  boiling  results  when  the 
phase  cancellation  is  not  complete;  the  residual  phase  encountered  with 
speckle  boiling  has  the  same  effect  on  the  correlation  function  as  phase 
aberration  has  on  the  impulse  response  function  of  an  imaging  system. 

2.3  CALCULATIONS  OF  NORMALIZED  CORRELATION 


In  this  section  we  apply  the  theory  developed  above  to  investigate 
speckle  from  rotating  objects.  For  simplicity  we  consider  the  2-D 
circular  disk  object  shown  in  Fig,  2-3.  Uniform  plane  wave  illumination 
is  used  with  the  direction  of  propagation  parallel  to  the  vector 
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P  x 


Figure  2-3.  Optical  system  used  in  the  simulation  to  investigate 
speckle  from  rotating  disk  objects. 


connecting  the  object's  axis  of  rotation  with  the  detection  point  p,. 
The  detection  point  p^  is  defined  to  be  located  at  the  p  =  0  position  on 
the  p-axis.  As  the  object  rotates,  the  speckle  pattern  in  the  detection 
plane  translates  from  p^  toward  p^.  For  the  calculations  presented  here 
we  took  the  object  radius  to  be  r  =  10  mm  and  the  range  z  =  100  mm. 

The  first  step  in  computing  the  speckle  correlation  function  is  to 
characterize  the  scattering  properties  of  the  object  by  specifying  the 
spatial  dependence  of  the  expected  scattering  cross-section  < I A ( x ) I  >. 
Here  we  consider  three  surface  types;  the  first  is  an  object  for  which 
the  scattering  cross-section  is  uniform  for  all  x.  This  would  result  if 
the  object  is  covered  with  perfect  retrorefl ecti ve  paint;  in  this  case 

<IA(x)12>u=1  ,  (2-11) 

where  the  subscript  U  stands  for  uniform.  Because  the  measured 
correlation  is  normalized,  Eq.  (2-11)  can  have  any  constant  value;  unity 
was  chosen  for  simplicity. 

The  second  type  of  surface  to  be  considered  is  Lambertian  f*.  r  which 
the  scattering  cross-section  varies  as  the  cosine  of  the  angle  between 
the  underlying  surface  normal  and  the  direction  of  observation.  For  the 
coordinate  system  shown  in  Fig.  2-3  with  z  »  r  we  have  approximately 

<IA(x)l2>L  =  [l  -  [f]2] 1/2  .  (2-12) 
where  the  subscript  L  stands  for  Lambertian  and  r  is  the  object  radius. 
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The  final  object  type  to  be  considered  has  a  Gaussian  cross-section 
function  given  by 

,  X21 

<  I A ( x )  r>f  =  exp  -  ,  (2-13) 

w 

where  the  subscript  G  stands  for  Gaussian.  For  the  simulations  reported 
here  with  a  10  mm  radius  object  we  used  w  =  5  mm  which  corresponds  to  an 
ooject  that  exhibits  largely  specular  reflection. 

The  correlation  functions  computed  using  Eq.  (2-10)  for  the  disk 
objects  are  shown  in  Figs.  2-4(a-c).  For  each  of  the  rotations 
considered  the  normalized  correlation  was  computed  for  several  locations 
of  P2;  thus  Fig.  2-4  shows  l^12'  as  a  function  of  detector  separation 
for  various  amounts  of  object  rotation,  A0  =  fl(t2  -  t^).  The 
correlation  for  each  point  on  the  curves  was  found  by  tracing  100  rays 
equally  spaced  along  the  x-axis.  The  wavelength  of  the  laser 
illumination  was  X  =  0.5  /*m. 

Figure  2-4(a)  contains  the  correlation  function  for  an  object  with 
uniform  scattering  cross-section  given  by  Eq.  (2-11).  Each  curve  shows 
only  the  most  significant  sidelobes.  The  first  curve,  A 6  =  0.0, 
designates  the  speckle  size;  the  first  zero  is  located  at  2.375  pm  while 
the  formula  Xz/2r  predicts  2.5  /im.  Notice  that  for  increased  rotation 
the  speckle  pattern  translates  and  boils. 

A  rough  estimate  for  the  amount  of  speckle  translation  is  well  known 
to  be  2zA0.  For  the  middle  curve  in  Fig.  2-4(a),  A0  =  0.2  mrad  and  the 
expected  translation  Is  thus  40.0  pm  while  the  observed  value  was  38.75 
^m.  For  A0  =  0.4  mrad  the  observed  translation  was  77.5  /*m  which  is 
precisely  twice  the  value  for  A0  =  0.2  mrad. 
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ia)  Uniform 


i  0.0  mrad 


-0  02  0  0.02  0.04  0.06  0.08  0.1 

Detector  Separation  (mm) 

(b)  Lambertian 
1  0  T  1  0.0  mrad 


-0.02  0  0.02  0.04  0.06  0  08  0.1 

Detector  Separation  (mm) 

(c)  Gaussian 

t  0.0  mrad  ,  n  0 


-0.02  0  0.02  0.04  0  06  0  03  0  1 

Detector  Separation  (mm) 

Figure  2-4.  Correlation  functions  for  objects  with  uniform, 

Lambertian,  and  Gaussian  scattering  profiles  are  shown  in 
(a),  (b)  and  (c)  respectively.  Each  plot  shows 
correlation  as  a  function  of  detector  separation  for  the 
various  amounts  of  object  rotation. 
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An  interesting  feature  of  Fig.  2-4(a)  is  the  sidelobes.  Notice  that 
their  magnitude  increases  steadily  with  rotation.  From  the  severity  of 
the  sidelobes  one  can  expect  that  the  phase  of  the  correlation  function 
given  in  Eq.  (2-10)  is  significantly  aberrated.  This  aberration  is 
further  studied  in  the  discussion  of  Fig.  2-5  below. 

Figure  2-4(b)  contains  the  correlation  function  calculated  for  the 
Lambertian  object  described  in  Eq.  (2-12).  The  same  trends  are  observed 
as  for  the  uniform  object,  however  now  the  speckle  size  is  slightly 
larger  and  the  boiling  and  sidelobes  are  less  severe. 

Figure  2-4(c)  shows  the  correlation  function  for  the  object  with 
Gaussian  cross-section  given  in  Eq.  (2-13)  with  parameter  w  =  5.0  mm. 
Notice  that  the  speckle  size  has  increased  over  the  previous  two 
surfaces  and  that  the  amount  of  boiling  and  sidelobe  magnitude  is 
decreased. 

To  gain  more  insight  Into  the  behavior  of  the  curves  shown  in 
Fig.  2-4,  we  have  plotted  the  phase  function  A^(x)  as  a  function  of  x  in 
Fig.  2-5  for  various  amounts  of  rotation.  For  all  cases  the  phase 
functions  plotted  are  for  the  separations  of  p1  and  p^  that  correspond 
to  the  correlation  peaks  in  Fig.  2-4(a).  We  have  also  subtracted  A^(0) 
and  thus  plotted  relative  phase.  Notice  that  the  aberration  increases 
with  the  amount  of  rotation.  To  understand  how  these  curves  affect  the 
speckle  dynamics  consider  Eq.  (2-10).  For  the  object  with  uniform 
cross-section  the  correlation  value  is  given  by  the  sum  over  the  phases 
from  the  entire  aberrated  wavefront  and  thus  with  increased  rotation  the 
speckle  boils  markedly  and  sidelobes  appear.  For  the  Lambertian  object 
the  contribution  from  the  highly  aberrated  edge  region  in  Eq.  (2-10)  is 
reduced  and  thus  the  correlation  is  higher  than  for  the  uniform  object. 
Finally  for  the  Gaussian  case  only  the  central  portion  of  the  wavefronts 
in  Fig.  2-5  contribute  and  thus  essentially  no  boiling  is  exhibited. 
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A  6  (x)  Waves 
0.5  r 


-10.0  -5.0  0  5.0  10.0 


Position  on  Object,  X  (mm) 


]n,r-a  - 5 .  The  phase  function  A^(x)  that  appears  in  Eq.  (2-5)  is 
plotted  as  a  function  of  position  on  the  object  for 
various  amounts  of  object  rotation. 
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One  important  observation  to  make  from  this  investigation  is  that 
the  surface  cross-section  has  the  same  effect  on  the  correlation 
function  as  an  apodizer  does  on  the  imoulse  response  of  an  imaging 
system. 

2.4  EXPERIMENTAL  MEASUREMENTS  OF  SPECKLE  CORRELATION 

In  this  section  results  of  experimental  measurement  of  the 
normalized  speckle  correlation  are  given.  The  optical  system  used  in 
the  experiments  is  shown  in  Fig.  2-6.  A  uniform,  collimated  and 
polarized  beam  from  an  Ar+  laser  operating  at  X  =  0.5145  /*m  is  reflected 
by  a  beamsplitter  and  illuminates  the  target.  The  target  was  a  4.6mm 
diameter  cylinder  with  height  4.6mm.  Three  different  surface  materials 
were  used;  they  are  summarized  in  Table  2.1.  These  materials  were 
chosen  to  demonstrate  that  the  theory  applies  to  a  broad  range  of 
materials.  A  linear  polarizer  was  placed  in  front  of  the  detector  to 
insure  that  a  single  polarization  was  recorded  even  If  the  target 
depolarized  the  illumination. 

The  speckle  patterns  were  recorded  using  a  CCD  camera  with  a  pixel 
size  of  18  pm  in  the  direction  of  speckle  translation  and  21.3  /*m  in  the 
other  direction.  The  camera  was  oriented  In  a  monostatic  arrangement  so 
that  the  central  pixel  was  on  the  virtual  axis  of  the  Illuminating  beam. 
The  distance  from  the  target's  axis  of  rotation  to  the  video  camera  was 
750  mm. 

The  output  from  the  video  camera  was  recorded  with  an  eight  bit 
digitizer.  Correlation  functions  were  computed  from  the  recorded 
speckle  patterns  using  an  array  processor.  The  relationship  between  the 
intensity  measurements  recorded  by  the  detector  array  and  the  normalized 
correlation  is  given  by  [2.9] 
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Figure  2-6.  Optical  system  used  in  the  experiments. 


22 


Material  A  Retro-reflective  paint  composed  of  glass 

beads  in  a  base  material.  3M  model  7216. 

Material  B  Typists  correction  fluid.  Produces  a 

uniform,  flat-white  surface. 

Material  C  Silver  paint. 


Table  2.1  Materials  used  to  coat  cylindrical  objects. 
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where  and  I2  denote  the  speckle  patterns  recorded  before  and  after 
rotation.  To  compute  the  ensemble  averages  in  Eq.  (2-14)  we  used  the 
statistically  stationary  nature  of  the  speckle  patterns  so  that  instead 
of  averaging  over  an  ensemble  of  rough  surfaces,  the  spatial  average  was 
computed.  Simulations  were  conducted  to  determine  that  the  speckle 
pattern  was  indeed  statistically  stationary  over  the  area  of  the 
recorded  speckle  patterns.  The  speckle  pattern  was  recorded  over  a  128 
x  128  pixel  array  both  before  and  after  object  rotation.  Each  of  these 
data  arrays  was  imbedded  in  a  256  x  256  array  of  zeros  and  the 
correlation  in  the  numerator  of  Eq.  (2-14)  was  computed  using  FFT 
techniques.  We  found  that  the  computed  values  of  |/*.2I  were  least 
noisy  if  the  expected  values  <^>  and  <I2>  in  Eq.  (2-14)  were  computed 
only  over  the  overlap  region  corresponding  to  each  correlation  offset. 
Thus  the  values  of  <I^>  and  <I2>  are  actually  arrays.  This  technique 
for  computing  l/<12l2  is  superior  to  using  constant  values  for  <^>  and 
<I2>  because  it  introduces  positive  correlation  between  the  numerator 
and  denominator  of  Eq.  (2-14)  which  enhances  the  signal -to-noise  ratio 
[2.11].  Each  of  the  correlation  functions  was  computed  from  a  single 
pair  of  speckle  data  frames;  there  was  no  averaging  over  multiple 
frames. 

Before  computing  the  normalized  correlation  we  had  to  correct  the 
speckle  data  for  pixel-to-pixel  non-uniformities  originating  at  the  CCD 
camera.  The  two  effects  we  considered  were  dark  current  and 
responslvity.  To  correct  the  responsivlty  we  recorded  a  frame  of  data 
with  the  sensor  exposed  to  uniform  light  from  a  distant  broadband 
incoherent  source  and  divided  each  frame  of  speckle  data  by  this 
correction  frame.  We  found  that  non-uniformities  in  the  dark  current 
were  insignificant,  however,  the  dark  current  was  responsible  for  adding 
a  uniform  bias  to  the  data.  To  correct  this  we  modified  the  bias  of  the 


24 


speckle  data  so  that  the  speckle  contrast  (the  ratio  of  the  standard 
deviation  to  the  mean)  was  unity.  Unity  contrast  is  a  property  of 
fully-developed  speckle,  so  by  appropriately  biasing  the  speckle  data 
the  correspondence  to  fully-developed  speckle  is  increased. 

Examples  of  experimentally  measured  correlation  functions  are  shown 
in  Figs.  2-7 (a)  and  (b) .  For  both  figures  the  data  was  collected  using 
the  setup  shown  in  rig.  2-6  and  the  object  was  coated  with  material  B. 
Figure  2-7(a)  corresponds  to  zero  object  rotation  and  in  Fig.  2-7 (b)  the 
object  has  rotated  2.0  arcmin.  These  photographs  result  from  recording 
the  128  x  128  frames  of  speckle  data  before  and  after  object  rotation 
and  computing  the  normalized  correlation  given  by  Eq.  (2-14).  The  grey 
levels  were  set  so  that  l/a^l  =  0.0  is  darkest  and  l/^l  =  1.0  is 
brightest  and  for  intermediate  values  the  grey  level  increases  linearly. 
If  the  measured  value  of  I n^\  was  below  zero  the  grey  level  was  set  to 
the  darkest  value  and  if  it  was  above  unity  the  grey  level  was  set  to 
the  brightest  value. 

The  central  pixel  in  Figs.  2-7 (a)  and  (b)  corresponds  to  zero 
correlation  offset;  thus  for  zero  rotation  (Fig.  2-7(a))  the  correlation 
peak  appears  in  the  center.  In  Fig.  2-7 (b)  the  object  has  rotated  2.0 
arcmin  and  the  corresponding  translation  of  the  speckle  pattern  gives 
the  shifted  correlation  peak. 

The  extent  of  the  correlation  functions  is  limited  by  the  128  x  128 
pixel  size  of  the  speckle  data  arrays.  Notice  that  the  statistical 
noise,  which  has  a  cloudlike  appearance,  increases  radially  from  the 
center  of  each  photograph.  This  is  because  the  number  of  independent 
speckles  used  to  compute  the  normalized  correlation  decreases  with 
increased  detector  offset.  As  a  result,  one  can  expect  false 
correlation  peaks  to  occur  in  the  corners  of  the  correlation  function 
where  relatively  few  speckles  enter  the  correlation  measurements. 
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EXPERIMENTAL  CORRELATION  MEASUREMENTS 
a)  0.0  arcmin  rotation  b)  2.0  arcmin  rotation 
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Figure  2-7.  Examples  of  the  speckle  correlation  functions  computed 

from  experimental  data,  (a)  Object  rotation  is  zero,  (b) 
Object  rotation  is  2.0  arc  min. 
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To  theoretically  predict  the  correlation  functions  for  the  various 
objects  we  used  the  ray  trace  procedure  discussed  in  Section  2.3  of  this 
paper.  Because  the  height  of  the  objects  was  small  (4.6mm),  we  found 
that  the  correlation  function  calculated  using  many  cross-sections  of 
the  object  was  identical,  in  the  direction  of  speckle  translation,  to 
the  correlation  function  computed  using  only  the  central  cross-section 
of  the  object.  In  fact,  we  found  that  the  height  of  the  cylinder  could 
be  increased  by  several  times  before  affecting  the  correlation  function. 
All  of  the  theoretical  correlation  functions  given  here  are  thus 
calculated  from  only  the  central  cross-section  of  the  object. 

To  compute  the  theoretical  correlation  function  we  had  to 
experimentally  measure  the  expected  scattering  profile,  < I A (x) 1  >,  for 
each  target  material.  To  accomplish  this  we  imaged  the  target  onto  the 
CCD  camera  by  placing  a  lens  between  the  CCD  camera  and  beamsplitter  and 
illuminated  the  target  with  spatially  incoherent  light  produced  by 
passing  the  laser  beam  through  a  rotating  diffuser.  This  procedure 
worked  well  for  the  objects  covered  with  materials  B  and  C,  however,  for 
material  A  more  averaging  was  needed  because  of  the  strong  point-like 
returns  from  the  glass  beads.  To  achieve  more  averaging,  we  rotated  the 
target  while  collecting  the  profile  data.  The  scattering  profiles  for 
target  materials  A,  B  and  C  are  shown  in  Fig.  2-8(a),  (b)  and  (c) 
respectively.  Notice  that  the  profiles  are  all  different  with  material 
A  being  the  most  diffuse  and  material  C  the  most  specular.  We  found 
that  the  computed  correlation  functions  are  dependent  on  gross  features 
of  the  profiles  and  insensitive  to  the  noise  in  Figs.  2-8(a-c). 

The  theoretically  calculated  correlation  functions  along  with  the 
experimental  measurements  are  shown  for  materials  A,  B  and  C  in  Figs.  2- 

9(a),  (b)  and  (c)  respectively.  The  correlation  functions  for  the 

2 

differjnt  materials  are  shown  in  separate  illustrations  with  I p^\ 
plotted  vs.  detector  separation.  The  solid  lines  are  theory  and  the 
circles  are  experimental  measurements.  The  experimental  results  are 
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cuts  in  the  direct: ~.i  of  speckle  translation  through  the  correlation 
peak  of  the  2-D  correlation  function  computed  by  the  array  processor  as 
shown  in  Figs.  2-7(a)  and  (b) .  In  each  figure,  the  correlation 
functions  are  shown  for  increasing  amounts  of  object  rotation  from  0  to 
3  arcmin.  Notice  that  for  zero  object  rotation  theory  and  experiment 
are  in  excellent  agreement  for  materials  A  and  B  while  the  agreement  for 
material  C  is  not  as  close.  For  material  C  the  experimental  correlation 
function  is  considerably  narrower  than  theory.  We  believe  that  this 
discrepancy  results  because  the  measured  scattering  profile  shown  in 
Fig.  2-8(c)  is  broader  than  the  true  profile  should  be.  This  indicates 
that  the  size  of  the  simulated  incoherent  source  (Illuminated  region  of 
rotating  diffuser)  was  too  large.  In  future  experiments,  especially  for 
objects  with  strong  specular  returns,  the  size  of  the  Incoherent  source 
should  be  better  controlled. 

With  increasing  object  rotation  the  speckle  translates  which  causes 
the  correlation  peaks  to  shift.  Theory  and  experiment  remain  in  good 
agreement  in  Figs.  2-9(a-c)  except  for  the  largest  rotation  for  which 
the  agreement  is  less  marked.  For  the  larger  rotations,  the  speckle  has 
translated  a  considerable  fraction  of  the  width  of  data  array  which 
gives  a  lower  signal-to-noise  ratio  that  is  evident  in  the  disagreement 
between  theory  and  experiment  regarding  the  peak  value  of  the 
correlation  function.  For  the  larger  rotations  also  notice  the  slight 
disagreement  between  theory  and  experiment  regarding  the  position  of  the 
correlation  peak.  At  this  point  we  do  not  know  where  the  disagreement 
originates.  It  seems  plausible  that  it  originates  in  measuring  the 
target  distance  and  rotation  or  the  detector  pixel  spacing,  but  we  were 
extremely  careful  in  measuring  these  parameters.  One  other  possible 
source  is  the  theory;  the  model  used  to  theoretically  compute  the 
correlation  functions  is  based  on  several  simplifying  approximations 
which  may  combine  to  cause  the  slight  discrepancy. 
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Figure  2-8.  Scattering  profiles  for  materials  A,  B  and  C  of  Table  1 
are  shown  in  (a),  (b)  and  (c)  respectively. 
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Figure  2-9.  Theoretical  and  experimental  correlation  functions  are 
compared  for  materials  A,  B  and  C  in  (a),  (b)  and  (c) 
respectively.  The  correlation  functions  are  plotted 
versus  detector  separation  for  the  various  amounts  of 
object  rotation.  Solid  lines  are  theory  and  circles  are 
experimental  measurements. 
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2.5  MULTIPLE  OBJECTS  ROTATING  AT  DIFFERENT  RATES 

Consider  the  dynamics  of  the  speckle  pattern  produced  when  multiple 
rotating  targets  are  illuminated  simultaneously.  Laboratory  observation 
of  such  a  speckle  pattern  reveals  that  each  rotation  rate  imparts  a 
component  of  translation  to  the  speckle  pattern,  however,  there  is 
considerable  boiling  because  the  patterns  from  each  target  add 
coherently. 

To  further  analyze  the  speckle  dynamics  consider  the  two  target 
case.  Eq.  (2-10)  gives 

H  <IA(xk)l2>  exp[iA^(xk)]  +  <IA‘(xk)l2>  exp[iA4' (xk)] 

„  =  J< - — - p - = - 2 -  (2'15) 

12  Z  < I A ( x  k ) I  /  +  H  <IA'(xk)l2> 

k  K  k  K 

where  the  primed  and  unprimed  values  of  A  and  denote  contributions 
from  the  two  targets.  In  Section  2.2  it  was  noted  that  for  a  single 
object  pure  speckle  translation  results  when  A^(xk)  is  constant  for  all 
xk,  and  the  amount  of  decorrelation  is  dependent  on  the  deviation  of 
A^(xk)  from  being  constant.  Based  on  this  observation  we  can  make  the 
following  conclusions  for  the  two  target  case: 

a)  if  A^(xk)  is  constant  for  both  targets  the  value  of  is  unity. 
This  situation  would  result,  for  example,  if  two  identical  targets 
were  at  the  same  distance  with  the  same  rotation  rate.  It  could 
also  result  if  the  two  targets  were  at  different  distances  with 
proportionately  different  rotation  rates. 

b)  If  A^(xk)  is  constant  for  one  object  and  deviates  markedly  for  the 
other  the  correlation  value  would  be 
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For  such  targets  we  can  conclude  that  a  separate  correlation  peak 
will  exist  for  each  and  the  corresponding  peak  value  is  given  by  the 
fraction  of  the  total  return  that  is  due  to  the  particular  object  as 
dictated  by  Eq.  (2-16).  Thus  for  two  targets  of  equal  intensity  the 
peak  values  are  ~  0.5. 

c)  If  the  values  of  A^x^)  are  neither  constant  nor  vary  substantially, 
Eq,  (2-15)  indicates  that  the  contributions  to  from  the  two 
objects  will  add  coherently. 

The  above  conclusions  also  extend  to  more  than  two  objects.  With 
several  objects  we  could  thus  expect  a  peak  for  each  distinct  rotation 
rate  with  the  value  at  each  peak  given  by  the  fraction  of  the  total 
return  from  the  specific  target. 

Experiments  were  performed  to  test  the  above  conclusions.  The 
experiments  were  identical  to  those  of  Section  2.4  except  that  two 
targets  were  used,  one  stationary  and  one  rotating.  The  Individual 
targets  were  identical  to  those  used  In  the  experiments  reported  in 
Section  2.4  with  cylinder  diameter  and  height  equal  to  4.6  mm.  Both 
were  coated  with  material  B.  The  rotating  object  was  placed  on  the 
optical  axis  of  the  system.  The  stationary  object  was  placed  above  the 
other  with  the  axes  of  symmetry  coaligned.  As  in  Section  2.4  we  found 
that  the  theoretical  predictions  based  on  a  single  cross  section  from 
each  object  were  equivalent  to  the  predictions  from  multiple  cross 
sections  and  therefore  single  cross-sections  were  used  to  form  the 
theoretical  predictions. 
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Theory  and  experiment  are  compared  in  Figs.  2-10(a-d).  In  each 

O 

figure  is  plotted  as  a  function  of  detector  separation  for 
rotation  values  of  0  arcmin  to  1.5  arcmin  in  increments  of  0.5  arcmin. 
The  solid  lines  are  theory  and  the  connected  circles  are  experimental 
measurements.  For  zero  rotation,  shown  in  Fia.  2-10(a),  theory  and 
experiment  are  in  excellent  agreement.  For  0.5  arcmin  of  rotation  in 
Fig.  2-10 (h )  there  is  excellent  agreement  regarding  the  position  of  the 
correlation  peaks,  however,  there  is  substantial  disagreement  regarding 
the  peak  values.  In  the  experiments  we  did  adjust  the  illumination  so 
that  the  peak  values  were  equal,  however,  the  peak  values  were  both 
roughly  =  0.40  in  Fig.  2- 10 (b)  while  the  theoretical  prediction 
is  slightly  greater  than  0.25  (for  larger  rotations  the  peak  values 
approach  a  minimum  of  0.25  as  discussed  above). 

For  larger  rotations  the  trends  continue;  excellent  agreement 
regarding  peak  location  with  experimental  peak  values  significantly 
larger  than  theory.  Also  notice  that  a  sidelobe  not  predicted  by  theory 
appears  in  Fig.  2-10(c).  In  Fig.  2~10(d)  the  correlation  peaks  are 
completely  separated. 

At  this  point  we  believe  that  the  experimental  values  are 
consistently  higher  than  theory  because  of  distortions  to  the  statistics 
of  the  speckle  patterns  imparted  by  the  CCD  camera  and  digitizer 
electronics.  In  the  experiments  the  contrast  of  the  speckle  patterns 
was  set  to  unity  which  Increases  the  correspondence  to  fully  developed 
speckle,  however,  we  believe  that  the  probability  density  functions  of 
the  data  deviated  from  being  negative  exponential.  In  future  research 
we  will  modify  the  instrumentation  in  hopes  of  improving  the  agreement 
between  theory  and  experiment. 
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;rrelation  functions  for  speckle  from  two  objects.  One 
Aject  is  stationary  and  the  other  object  rotates 
increasing  amounts  in  (a),  (b),  (c)  and  (d)  as  shown. 

;  1  id  lines  are  theory  and  connected  circles  are 
experimental  measurements. 
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2.6  COMPLICATED  ROTATING  OBJECTS 

Experiments  were  also  conducted  using  objects  with  complicated 
underlying  shape.  The  experimental  setup  shown  in  Fig.  2-6  was  used  and 
the  object  was  a  metallic  scale  model  of  a  military  jeep.  The  surfaces 
of  the  object  were  rough  and  the  underlying  shape  had  several  concave 
features.  The  dimensions  of  the  model  in  the  directions  transverse  to 
the  illumination  were  roughly  15  mm  by  8  mm  so  the  speckle  size  was 
considerably  smaller  than  for  the  cylinders  used  in  Section  2-4.  The 
results  of  the  experiments  exhibited  the  same  trends  as  for  the 
cylindrical  objects;  the  speckle  translation  is  proportional  to  object 
rotation  with  negligible  boiling.  The  only  detectable  difference  wes 
the  speckle  size  which  resulted  in  correlation  peaks  that  were 
considerably  narrower.  From  these  results  one  can  conclude  that 
rotation  measurements  can  be  accomplished  equally  well  for  objects  with 
simple  or  reasonably  complicated  underlying  shapes. 

2.7  ROTATION  MEASUREMENT  AT  LOW  LIGHT  LEVEL 

In  this  section  we  analyze  the  ability  to  perform  rotation 
measurement  at  low-light  level.  The  SNR  for  performing  low-light  level 
speckle  correlation  has  been  considered  extensively  in  the  literature. 
Analysis  applicable  to  this  application  is  contained  in  Ref.  2.12.  Let 
us  adapt  the  results  of  Ref.  2.12  to  consider  the  case  using  a  2-D 
detector  array  to  measure  rotation.  Let  us  assume  that  the  speckle 
pattern  does  not  move  in  the  tir.e  between  sequential  recordings  of  the 
2-D  speckle  pattern.  As  shown  in  previous  sections,  we  are  thus 
concerned  with  measuring  a  correlation  function  that  has  a  peak  value  of 
l/il  =  1.0  at  the  origin  and  gradually  goes  to  zero  value  in  a  region 
equal  to  the  speckle  size  of  XR/D.  The  SNR  for  such  correlation 
measurements  is  [2.12] 
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4MK 


il/2 


20  +  12/N<n>  +  2/N<n/J 


(2-17) 


where  M  is  the  number  of  independent  speckles  encompassed  by  the 
detector  array,  K  is  the  number  of  independent  frame  pairs  used  to 
measure  the  correlation,  N  is  the  number  of  detection  channels  per 
speckle  and  <n>  is  the  mean  number  of  photons  per  detection  channel. 
With  the  number  of  photons  per  speckle  large  (N<n>  »  1),  it  follows 
that  the  first  term  in  the  denominator  dominates.  This  first  term  is 
thus  associated  with  statistical  noise  imparted  by  the  speckle.  As  the 
number  of  photons  per  speckle  decreases,  the  last  two  terms  in  the 
denominator  of  Eq.  (2-17)  dominate.  To  minimize  noise  at  these  low 
light  levels,  it  follows  by  inspection  of  the  last  terms  that  with  N<n> 
being  constant,  the  number  of  detectors  per  speckle  N,  should  be  as 
small  as  possible. 


We  conducted  some  low-light  level  simulations  to  test  the  validity 
of  the  SNR  formula  in  Eq.  (2-17).  In  the  laboratory  we  recorded  a  high¬ 
light  level  speckle  pattern  produced  by  back-illuminating  a  small 
diffuse  square  aperture  with  laser  light.  The  pattern  was  recorded 
using  a  2-D  detector  array.  A  region  of  size  128  x  128  pixels  was  used 
in  the  simulation.  The  detector  pixel  spacing  and  illuminating  aperture 
size  were  such  that  there  were  9  and  11  pixels  per  speckle  in  the 
vertical  and  horizontal  directions  respectively.  These  dimensions  were 
extracted  by  correlating  the  high  light  level  speckle  pattern.  In  this 
simulation  we  thus  had  M  =  164  and  N  =  99. 


To  simulate  speckle  recorded  at  low  light  level,  a  Poisson  random 
number  generator  was  applied  to  the  recorded  speckle  pattern.  For  the 
Poisson  generator  we  used,  the  expected  number  of  photons  per  pixel, 
<n>,  is  equal  to  the  grey  value  per  pixel  averaged  over  the  recorded 
pattern.  Thus  to  achieve  specific  values  of  <n>  the  recorded  pattern 
was  multiplied  by  a  scaling  constant  and  then  passed  through  the  Poisson 
random  number  generator. 
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One  goal  of  these  simulations  is  to  determine  the  lowest  light  level 
at  which  the  rotation  measurements  can  be  made.  In  the  limit  of  low 
light  level,  the  expression  for  the  SNR  reduces  to 


SNR 


(2-18) 


where  we  have  taken  K  =  1.  By  setting  the  SNR  to  unity  and 

substituting  M  =  164  and  N  =  99,  we  find  that  the  lower  limit  of 

-2 

acceptable  light  level  is  <n>  =  5.54  x  10  photons  per  pixel,  or  0.55 
photons  per  speckle.  Note  that  if  the  experiment  conducted  with  1 
pixel  per  speckle  (N  =  1),  the  lowest  light  level  would  be  0.055 
photons  per  speckle.  Thus  the  lowest  acceptable  light  level  is 
decreased  by  an  order  of  magnitude  by  not  oversampling  this  speckle. 
Oversampling  should  therefore  be  avoided  in  experimental  design. 

The  results  of  this  simulation  are  shown  in  Figs.  2-11,  2-12  a  ; 
2-13.  Figure  2-11  shows  the  original  speckle  pattern  and  the  pattern 
at  levels  of  104 ,  103,  10^,  10,  5,  3  and  1  photons  per  speckle.  Notice 
that  the  appearance  of  the  speckle  pattern  gradually  degrades  as  the 
light  level  decreases. 

To  simulate  the  correlation  measurements  at  low  light  level  we 
generated  two  independent  realizations  of  the  speckle  noise  for  both 
from  the  same  digitized  speckle  pattern.  These  were  then  correlated 
using  the  same  Fourier  transform  procedure  as  in  the  experiments 
reported  in  Section  2.4.  The  simulation  thus  corresponds  to  situations 
in  which  the  object  has  not  rotated  between  the  recordings  of  the 
speckle  pattern. 
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SPECKLE  PATTERNS  FOR  URRIOUS  LIGHT  LBJELS  (PHOTOHS/SPEDCLO 


'  'jure  2-11.  Examples  of  speckle  patterns  are  shown  for  the  jr 
light  levels. 


Figure  2-12.  Correlation  functions  genera* 
1 evel s . 
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oVjts  taken  horizontally  through  the 

-■  .  •-•-•v  V- :  on  functions  shown  in  Fig.  2-12. 
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The  correlation  functions  obtained  for  the  various  light  levels  are 
shown  in  Fig.  2-12.  Along  with  the  high  light  level  correlation 
function  the  figure  shows  the  correlation  function  at  light  levels  of 
104,  103 ,  102,  10,  5,  3  and  1.  Note  that  the  correlation  peak  is 
recognizable  at  3  photons  per  speckle  and  that  it  is  not  recognizable 
at  1  photon  per  speckle. 

To  better  characterize  the  correlation  performance,  horizontal 
cross-sections  through  the  correlation  peak  are  shown  for  the  various 
light  levels  in  Figs.  2-13(a)-(g).  Even  at  100  photons  per  speckle  in 
Fig.  2- 13 (c) ,  the  correlation  function  is  barely  degraded.  The 
degradation  increases,  however,  until  it  is  not  distinguishable  in  Fig. 
2-13 (g)  at  1  photon  per  speckle. 

For  these  simulations  we  have  found  that  the  lowest  acceptable 
light  level  is  roughly  1  photon  per  speckle.  This  result  Is  in  general 
agreement  with  the  result  of  0.5  photons  per  speckle  predicted  by  Eq. 
(2-18).  It  must  be  stressed  that  this  lower  limit  on  the  light  level 
is  for  specific  values  of  M  and  N  and  for  other  configurations  the 
lower  limit  on  the  light  level  can  vary  substantially. 

2.8  COMPARISON  TO  CONVENTIONAL  DOPPLER  RADAR 

The  basic  difference  between  rotation  rate  measurement  using 
speckle  intensity  and  conventional  Doppler  radar  is  that  the  former 
uses  autodyne  detection  while  the  later  uses  heterodyne  detection. 
With  autodyne  detection  the  signal  results  from  the  interference  of 
light  scattered  from  different  regions  of  the  object  and  the  phenomenon 
of  Doppler  beating  corresponds  to  speckles  passing  over  the  detector. 
For  heterodyne  detection  the  signal  corresponds  to  the  interference  of 
light  from  the  entire  object  with  a  local  oscillator  and  the  Doppler 
signal  is  caused  primarily  by  relative  motion  between  the  transmit/ 
receive  platform  and  the  object.  The  heterodyne  Doppler  signal  is  also 
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modulated  by  the  speckles  passing  over  the  detector.  In  fact,  this 
speckle  modulation,  which  gives  rise  to  the  Doppler  spread,  is  the 
signal  component  that  gives  the  rotation  rate  information. 

In  a  comparison  conducted  in  Ref.  2.8,  it  was  found  that  autodyning 
and  heterodyning  are  equally  sensitive  at  measuring  Doppler  spread  and 
thus  rotation  rate.  It  was  also  shown  that  autodyning  is  far  less 
sensitive  to  the  optical  quality  of  the  receiver  and  laser  frequency 
stability.  For  SOI  applications  laser  frequency  stability  over  the 
long  propagation  distances  is  a  fundamental  problem  for  heterodyne 
systems. 

One  other  issue  is  the  SNR  for  autodyne  detection  vs.  that  for 
heterodyne  detection.  Let  us  assume  that  the  return  signal  is  shot 
noise  limited;  this  assumption  is  applicable  to  the  SDI  scenario. 
Noise  sources  such  as  background  noise  and  dark  current  are  thus  taken 
to  be  small  compared  to  the  signal  level.  It  follows  that  In  this 
limit  the  SNR's  for  the  two  detection  schemes  are  equivalent;  both  are 
capable  of  single  photon  detection  [2.13].  If  the  return  signal  is  not 
shot  noise  limited,  heterodyne  detection  does  allow  one  to  achieve  the 
shot  noise  1 imit  [2.13] . 

In  summary,  speckle  intensity  methods  and  heterodyne  methods  can, 
in  theory  perform  equally  well  in  measuring  object  rotation  rate.  In 
practice,  however,  the  heterodyne  methods  require  optical 
instrumentation  of  much  higher  quality  and  lasers  with  much  greater 
frequency  stability.  These  requirements  are  particularly  imposing  when 
a  large  array  of  heterodyne  receivers  is  required.  Using  current 
technology,  we  do  not  feel  that  heterodyne  detection  is  feasible  for 
the  SDI  scenario,  however,  it  does  provide  valuable  information, 
particularly  for  image  reconstruction,  and  does  warrant  technology 
development. 
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MEASUREMENT  OF  OBJECT  SEPARATION  USING  LASER  SPECKLE 


3.1  INTRODUCTION 

Consider  a  distant  object  consisting  of  two  separate  segments  as 
shown  in  Fig.  3-1 (a).  The  ability  to  resolve  the  segments  using  a 
conventional  imaging  system  can  be  characterized  by  the  Rayleigh 
criterion  and  the  minimum  resolvable  angular  separation  is 

M  =  1.22  X/D  ,  (3-1) 

where  X  is  the  wavelength  and  D  is  the  diameter  of  the  imaging  lens.  It 
follows  that  to  measure  the  separation,  rather  than  simply  resolve  it, 
requires  a  significantly  larger  aperture  [3.1].  The  achievable 
measurement  accuracy  is  thus  limited  by  the  ability  to  construct  large 
aperture,  diffraction  limited  optics. 

In  this  section  of  this  report  we  analyze  an  alternative  method  for 
separation  measurement  that  is  based  on  laser  speckle.  This  method 
requires  active  laser  illumination  of  the  separated  object  which  is 
assumed  to  be  diffusely  reflective.  Rather  than  using  an  imaging  lens, 
the  speckle  pattern  is  detected  using  a  2-D  detector  array.  While  this 
technique  does  not  relax  the  size  requirement  of  the  detection  aperture, 
it  does  eliminate  the  need  for  large  aperture,  diffraction  limited 
optics  by  instead  using  a  2-D  detector  array. 

The  basic  idea  behind  this  speckle  technique  is  illustrated  in  Fig. 
3-1.  An  example  of  a  speckle  pattern  from  the  separated  object  shown  in 
Fig.  3-1  (a )  is  shown  in  (b).  Notice  the  spatial  modulation  of  the 
speckle  pattern  introduced  by  the  object  separation.  The  frequency  of 
this  modulation  is  proportional  to  the  amount  of  separation.  This 
spatial  frequency  information  can  be  extracted  using  the  Fourier 
transform.  The  squared  magnitude  of  the  Fourier  transform  of  the 
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'jure  3-1.  Separation  measurement  using  laser  speckle,  (a)  Separated 
object,  (b)  Realization  of  speckle  pattern.  (c)  Fourier 
transform  of  speckle  intensity. 
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speckle  intensity  is  shown  in  Fig.  3-1 (c) ;  this  transform  essentially 
gives  a  speckled  version  of  the  object's  autocorrelation.  Measurement 
of  the  object  separation  is  reduced  to  computing  the  distance  between 
the  center  of  this  transform  and  the  adjacent  lobes.  In  Section  3.3  we 
discuss  the  use  of  centroid  measurements  to  extract  this  information. 

In  summary,  this  investigation  concerns  a  technique  for  measurement 
of  target  separation  using  laser  speckle.  While  this  technique  does  not 
reduce  the  aperture  requirement  below  that  of  a  conventional  imaging 
system,  it  does  not  require  a  diffraction  limited  lens  to  cover  the 
aperture,  but  instead  a  detector  array.  The  requirement  for  precision 
optics  is  thus  reduced. 

In  Section  3.2  of  this  report,  statistical  properties  of  the  Fourier 
transform  of  a  speckle  pattern  are  investigated.  The  relationship 
between  the  transformed  speckle  pattern  and  the  autocorrelation  of  the 
incoherent  representation  of  the  object  is  derived.  A  brief  discussion 
of  the  sampling  frequency  required  to  record  the  speckle  pattern  without 
aliasing  is  also  given. 

In  Section  3.3  the  measurement  of  object  separation  by  computing 
centroids  of  the  transformed  speckle  pattern  is  examined.  The 
measurement  accuracy  is  quantified  using  laboratory  data  by  comparing 
the  actual  and  measured  separations  for  several  realizations  of  the 
speckle  pattern  with  varying  numbers  of  speckles  over  the  detection 
aperture. 

3.2  PROPERTIES  OF  THE  FOURIER  TRANSFORM  OF  A  SPECKLE  PATTERN 

Let  I(x)  be  the  intensity  of  a  speckle  pattern  produced  by 
illuminating  a  distant,  rough  object  with  coherent  light.  To  simplify 
notation,  this  treatment  is  conducted  in  one  dimension,  however,  the 
analysis  is  directly  extendable  to  the  2-D  case.  It  is  assumed  that 
I (x)  represents  the  intensity  of  a  fully-developed  speckle  pattern  in 
which  case  I (x)  is  a  negative  exponential  random  variable. 
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■*'  transform  of  an  ideally  recorded  version  of  the  speckle 

liat*  1  "ten  as 

b/2 

I(u)  =  j  I (x)  exp(-i2r  ux)  dx  ,  (3-2) 

-b/2 

where  b  denotes  the  extent  of  the  detection  aperture.  Inspection  of  Eq. 
(3-2)  reveals  that  for  each  value  of  u,  with  the  exception  of  u  near 
zero,  ?  is  essentially  a  sum  of  randomly  weighted  phasors.  Excluding  u 
near  zero,  it  follows  from  central  limit  theorem  that  when  the  detector 
array  encompasses  many  independent  speckles,  I  is  a  circular,  complex 
Gaussian  random  variable.  The  Fourier  transform  of  the  speckle 
intensity  thus  obeys  the  same  statistical  model  as  the  optical  field  in 
a  speckle  pattern. 

The  expected  value  of  T  is  written  as 

b/2 

<T(u)>  =  f  <I(x)>  exp(-i 2r  ux)  dx  ,  (3-3) 

-b/2 

where  the  angular  brackets  denote  an  ensemble  average.  Because  the 
object  is  taken  to  be  at  a  large  distance  from  the  detection  plane,  it 
can  readily  be  assumed  that  the  speckle  pattern  is  spatially  stationary, 
in  which  case  <I(x)>  is  independent  of  x.  By  defining  <I(x)>  =  <I>  we 
have 

b/2 

<I(u)>  =  <I>  J  exp(-i2ir  ux)  dx  , 

-b/2 

=  b<I>  sinc(bu)  ,  (3-4) 


48 


JtERIM 


where  sinc(x)  =  (sin  tx)/tx.  In  the  limit  of  large  b  it  follows  that 

Lim  <I(u)>  =  <I>  5 (u)  ,  (3-5) 

b+® 


where  <5  is  the  Dirac  delta  function.  Thus  the  transform  of  the 
intensity  has  zero  expected  value  except  near  u  =  0. 

The  expected  value  of  the  squared  magnitude  of  1,  or  power  spectrum, 
is  written  as 


b/2  b/2 

<lT(u)l2>  =  f  I  <I(x1)  I(x2)>  exp[-i2ru(x1  -  x?)]  dXj  dx2  (3-6) 
-b/2  -b/2 

The  expected  value  that  appears  on  the  RHS  of  Eq.  (3-6)  can  be 
factorized  as  follows  [3.2] 

<I(x1)  I(x2)>  =  <I>2  -  x2)l2  +  1)  ,  (3-7) 

where  n  is  the  complex  normalized  correlation  of  the  speckle  pattern 
which  has  a  value  in  the  range  0  i  l/*l  $  1.  As  noted  in  Ref.  3.2,  p  is 
given  by  the  normalized  Fourier  tran'  'orm  of  the  intensity  distribution 
of  the  object.  This  relationship  is  analogous  to  the  van  Cittert- 
Zernike  theorem  from  classical  coherence  theory  [3.3].  In  writing  Eq. 
(3-7)  we  have  used  the  stationary  nature  of  the  speckle  pattern  and  thus 
only  the  difference  of  spatial  coordinates  appears. 

To  further  evaluate  Eq.  (3-6)  with  Eq.  (3-7)  substituted,  one  must 
transform  the  integration  to  the  difference  coordinates,  r  =  x^  -  x2, 
which  gives 
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b 

<  1 1  (u)  1 2>  =  <I>2  J  (b  -  Irl)  (l^(r)l2  +  1)  exp(-i2ir  ur)  dr  .  (3-8) 

-b 

Notice  that  Eq.  (3-8)  is  the  Fourier  transform  of  a  product.  The 
important  components  of  this  product  are  first 

b 

j  (b  -  Irl)  exp(-i2r  ur)  dr  =  b2  sinc2(bu)  ,  (3-9) 

-b 

which  by  assuming  that  b  is  much  greater  than  the  speckle  size  gives 

b  * 

f  l/i(r)  1 2  exp(-12r  ur)  dr  =  J  f({)  f(f  -  u)  d£  ,  (3-10) 

-b  -* 

'  =  f(u)  *  f(u)  (3-11) 

where  f  is  the  Fourier  transform  of  p  and  *  denotes  autocorrelation.  As 
mentioned  above  it  follows  from  the  van  Ci ttert-Zerni ke  theorem  that  f 
is  the  intensity  distribution  of  the  object.  Using  Eqs.  (3-9),  (3-10) 
and  (3-11),  Eq.  (3-8)  is  re-written  as 

<II(u)l2>  =  b2  <I>2  sinc2(bu)  *  (.f(u)  *  f (u) ]  +  b2<I>2  sinc2(bu)  (3-12) 

The  most  important  term  in  Eq.  (3-12)  is  the  autocorrelation  of  the 
object's  intensity  distribution,  f  *  f.  As  shown  in  Fig.  3-1,  this 
autocorrelation  gives  the  distinct  lobes  from  which  the  separation  is 
measured. 

Added  to  the  autocorrelation  is  a  sine  function  which  for  large  b 

reduces  to  a  delta  function  at  the  origin.  The  result  of  the  addition 

2 

is  then  convolved  with  a  sine  that  also  amounts  to  a  delta  function  for 
larbe  b. 
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One  other  important  feature  of  the  Fourier  transform  of  the  speckle 
pattern  is  the  size  of  the  speckle  in  the  transform  domain.  If  the 
original  speckle  data  is  contained  in  an  array  of  length  M,  which  is 
then  imbedded  in  a  larger  array  of  zeros  of  length  N,  the  speckle  size 
in  the  transform  domain  is  N/M. 

One  final  consideration  is  how  finely  the  speckle  pattern  must  be 
sampled  to  avoid  aliasing  in  the  transform  of  the  speckle  intensity. 
Let  us  define  d  to  be  the  size  of  the  1-D  object;  that  is,  the  object 
intensity  is  zero  outside  of  a  region  of  size  d.  It  follows  that  the 
maximum  detector  pixel  spacing  before  aliasing  occurs  in  the  transform 
i  s 


s  =  XZ/2d  (3-13) 

where  Z  is  distance  between  the  object  and  detector  array. 

3.3  MEASUREMENT  OF  OBJECT  SEPARATION 

In  the  last  section  it  was  shown  that  the  Fourier  transform  of  the 
intensity  of  a  speckle  pattern  is  essentially  a  speckled  version  of  the 
autocorrelation  of  the  object's  intensity  distribution.  In  this  section 
we  consider  the  extraction  of  object  separation  from  the 
autocorrelation. 

One  way  to  measure  object  separating  from  the  autocorrelation  would 
be  to  measure  the  distance  from  the  peak  of  one  of  the  sidelobes  to  the 
center  (see  Fig.  3- 1 (c) ) .  Because  of  the  speckle,  however,  this 
procedure  is  likely  to  be  inaccurate  especially  for  objects  with 
smoothly  peaked  autocorrelation  sidelobes. 

Another  method  for  separation  measurement  is  to  compute  the  centroid 
of  the  sidelobe  which  is  defined  by 
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x  =  Z  m  I Iml2  ,  (3-14) 
m 

*nere  M  is  the  sidelobe  mass 

M  =  H  IImi2  .  (3-15) 
m 

The  sums  in  Eqs.  (3-14)  and  (3-15)  extend  only  over  the  area  of  the 
sidelobe.  Because  the  center  of  the  autocorrelation  is  at  m  =  0,  the 
centroid  gives  the  object  separation  directly.  Two  dimensional 
separations  are  computed  by  defining  the  sums  in  Eqs.  (3-14)  and  (3-15) 
to  extend  over  the  2-D  sidelobe  and  y  is  computed  analogously  to  x. 

For  objects  with  diameter  larger  than  their  separation  the 
sidelobes  are  not  separated.  In  this  case  one  can  use  a  thresholding 
procedure  whereby  the  threshold  is  raised  until  the  sidelobes  are 
distinct  and  then  the  centroid  is  computed. 

Analytic  computation  of  the  accuracy  for  separation  measurement 
using  centroiding  is  very  complicated.  For  this  reason  we  have 
determined  the  accuracy  empirically.  Particularly  we  have  investigated 
the  number  of  speckles  encompassed  by  the  detector  array  required  to 
accurately  measure  the  separation.  This  aperture  requirement  can  then 
be  compared  to  the  Rayleigh  criterion  given  in  Eq.  (3-1)  and  the 
analysis  of  Cef.  i. 


Experimental  separation  measurements  were  conducted  using  the  setup 
shown  in  Fig.  3-2.  The  diffuse  object  consisted  of  two  vertically 
separated  circles  both  of  diameter  3.05  mm  with  center-to-center  spacing 
equal  to  twice  the  diameter.  This  object  was  illuminated  with  coherent 
illumination  at  X  =  0.5145  ^m  and  the  speckle  pattern  was  recorded  at  a 
distance  of  154  cm  using  a  CCD  detector  array.  The  center-to-center 
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Figure  3-2.  Experimental  setup  used  for  separation  measurements. 
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pixel  spacing  of  the  detector  array  was  Sx  =  21.3  pm  in  the  horizontal 
direction  and  =  18.0  pm  in  the  vertical  direction.  These  pixel 
'Oacings  are  well  under  the  upper  limits  established  in  Eq.  (3-13;. 

If,  for  example,  the  autocorrelation  sidelobe  has  a  centroid  of  y  in 
the  vertical  direction,  the  angular  separation  of  the  objects  in  the 
vertical  direction  is 

N 

«„  ■  s .  (3-16) 
y  — y  Ly 

where  is  the  number  of  detector  elements  in  the  vertical  direction  of 
the  array  including  zero  padding.  For  example,  if  the  speckle  data 
array  has  64  x  64  pixels  and  is  imbedded  in  a  256  x  256  array,  the  value 
of  L,  is  256. 

y 

In  the  experiments,  we  took  sections  of  the  speckle  data  with  sizes 
of  128  x  128,  64  x  64,  and  32  x  32  pixels  and  imbedded  this  data  into  an 
array  of  zeros  with  size  256  x  256  pixels.  The  left  column  of  Fig.  3-3 
shows  examples  of  the  speckle  data  for  the  different  data  array  sizes. 
The  center  and  right  columns  show  the  autocorrelation  estimates  for  two 
independent  realizations  of  the  speckle  data.  Notice  that  the  quality 
of  the  autocorrelations  decrease  as  the  size  of  the  speckle  data  array 
decreases . 

The  autocorrelation  centroids  were  computed  for  the  various  speckle 
data  array  sizes  and  were  compared  to  the  expected  centroids.  The 
results  of  this  experiment  are  shown  in  Fig.  3-4.  From  Eq.  (3-16)  the 
expected  object  separation  is  y  =  35.4  pixels.  This  is  shown  as  the 
horizontal  line  in  the  Fig.  3-4.  The  horizontal  axis  of  Fig.  3-4  gives 
the  dimension  of  the  speckle  data  array  size.  For  each  of  the  array 
sizes  (32,  64  and  128),  the  centroid  was  computed  for  three  realizations 
of  the  speckle  data.  Notice  that  the  accuracy  increases  with  the  size 
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Figure  3-3.  Examples  of  various  sized  speckle  patterns  are  shown  in 

the  left  column.  Two  realizations  of  the  autocorrelation 
from  independent  speckle  patterns  are  shown  in  the  center 
and  right  columns. 
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Figure  3-4.  A  comparison  of  experimentally  measured  separation  to  the 
true  value  for  various  speckle  pattern  sizes.  For  each 
size  3  independent  realizations  of  the  separation 
measurement  are  shown. 


of  the  data  array.  For  the  largest  array  size,  all  of  the  separation 
measurements  appear  within  one  pixel  of  the  true  value.  For  comparison, 
the  smallest  array  size  (32)  is  roughly  4  times  the  size  of  the  minimum 
aperture  required  to  resolve  the  separation  as  dictated  by  the  Rayleigh 
criterion  in  Eq.  (3-1) . 
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PHASE  RETRIEVAL  FOR  A  COMPLEX-VALUED  OBJECT 
USING  A  LOW-RESOLUTION  IMAGE 


4.1  SUMMARY 

It  is  very  difficult  to  reconstruct  an  image  of  a  complex-valued 
object  from  the  modulus  of  its  Fourier  transform  (i.e.,  retrieve  the 
Fourier  phase)  except  for  some  special  cases.  Using  additionally  a 
low-resolution  intensity  image  from  a  telescope  with  a  small  aperture, 
a  fine-resolution  image  of  a  general  object  can  be  reconstructed  by  a 
two-step  approach.  First  the  Fourier  phase  over  the  small  aperture  is 
retrieved  using  the  Gerchberg-Saxton  algorithm.  Then  that  phase  is 
used,  in  conjunction  with  the  Fourier  modulus  data  over  a  large 
aperture  together  with  a  support  constraint  on  the  object,  to 
reconstruct  a  fine-resolution  image  (retrieve  the  phase  over  the  large 
aperture)  by  the  iterative  Fourier  transform  algorithm. 

4.2  BACKGROUND 

Phase  retrieval  from  a  single  intensity  distribution  for  the  case 
of  compl ex-val ueH  objects  arises  in  a  number  of  applications  such  as 
holography,  wavefront  sensing,  and  imaging  with  coherent  illumination. 
If  the  support  of  the  object  (the  set  of  points  over  which  it  is 
nonzero)  is  well  known  or  of  a  favorable  type,  then  it  is  often 
possible  to  reconstiuct  an  image  of  the  object  from  the  modulus  of  its 
Fourier  transform  (the  square  root  of  the  Fourier  intensity)  using  the 
iterative  Fourier  transform  algorithm  [4.1].  Favorable  support 
constraints  include  polygons  with  no  parallel  sides  (particularly 
triangles),  which  must  be  known  a  priori  [4. 1,4.2],  and  supports  with 
separated  parts,  which  need  not  be  known  a  priori  [4. 1,4.3],  If,  on 
the  other  hand,  the  object  hac  a  polygonal  support  that  is  rot  known 
very  well  a  priori,  or  if  the  object  has  tapered  edges,  then  both  the 
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ability  to  converge  to  a  solution  and  the  quality  of  the  reconstructed 
image  deteriorate  [4. 1 , 4 .4-4 . 6] .  This  contrasts  sharply  with  the  case 
of  real,  nonnegative  objects,  for  which  phase  retrieval  is  much  easier 
[4. 7-4. 9].  Image  reconstruct i on  is  also  possible  if  the  complex-val ued 
object  has  a  strong  glint  or  glints  (for  example,  a  single  glint  well- 
separated  from  the  object  gives  rise  to  a  hologram  which  can  be 
reconstructed  easily  [4.10]). 

In  this  chapter  we  show  that  even  the  difficult  types  of  complex¬ 
valued  objects  can  be  reconstructed  if  one  has  a  low-resolution 
intensity  image  of  the  object,  taken  through  a  telescope  having  a  small 
aperture  contiguous  with  the  Fourier  intensity  measurements,  to 
supplement  the  Fourier  intensity  data.  I n  Section  4.3  an  example  is 
given  of  an  optical  system  that  would  produce  the  desired  measurements. 
Section  4.4  describes  the  data  processing  steps  required  to  reconstruct 
a  fine-resolution  image.  A  two-step  method  is  used,  employing  an 
accelerated  version  of  the  Gerchberg-Saxton  algorithm  [4.11-4.15]  in 
the  first  step  and  a  modified  version  of  the  iterative  Fourier 
transform  algorithm  [4.1,4,7-4.9,4.13-4.15]  in  a  second  step.  This  new 
modification,  the  expanding  weighted  modulus  algorithm,  was  necessary 
for  convergence  with  a  reasonable  number  of  iterations.  In  Section  4.5 
an  example  of  reconstruct! ng  an  image  using  this  approach  is  given,  and 
in  Section  4.6  are  conclusions. 

4.3  OPTICAL  SENSOR  CONFIGURATION 

Suppose  that  the  object  being  imaged  is  illuminated  by  a  coherent 
laser  and  is  far  away  so  that  the  relationship  between  the  optical 
field  at  the  object,  f(x),  and  that  in  the  aperture  plane  of  the 
optical  receiver,  f(u),  is  approximately  a  Fourier  transform  [4.16]. 
Here  u  and  x  are  both  two-dimensional  coordinates:  u  in  the  aperture 
plane  and  /  in  the  object  or  image  plane.  (If  it  is  a  Fresnel 
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transform,  then  the  method  described  here  will  work  with  minor 
modifications  [4.17].)  Figure  4-1  depicts  an  example  of  an  optical 
receiver  that  gathers  the  type  of  data  needed  for  the  reconstruction 
described  here.  An  array  of  light-bucket  detectors  (shown  with  field 
lenses  in  front  of  them)  samples  the  intensity  of  the  optical  field  in 
the  aperture  plane.  In  order  to  adequately  sample  the  intensity  of  the 
speckle  pattern  in  the  aperture  plane,  there  must  be  at  least  two 
detectors  per  speckle  width  in  each  dimension  (as  determined  by  the 
wavelength  of  the  laser,  the  distance  to  the  object,  and  the  object 
diameter).  Since  only  the  intensity  is  detected,  these  measurements 
are  independent  of  any  phase  errors  that  may  be  present  due  to 
atmospheric  turbulence  in  front  of  the  aperture  (assuming  that 
atmospheric  scintillation  is  negligible)  or  due  to  misalignment  of  the 
detector  array.  In  addition,  imbedded  in  the  array  (or  contiguous  with 
the  array  on  the  edge  of  the  array)  is  a  smal 1 -aperture  diffraction- 
limited  telescope.  If  located  on  earth,  the  smal 1 -aperture  telescope 
could  be  diffraction-limited  by  virtue  of  having  an  adaptive  optics 
system  that  compensates  for  atmospheric  turbulence  in  real  time.  Such 
an  adaptive  optics  system  may  not  be  practical  for  a  telescope  with  an 
aperture  the  size  of  the  entire  large  aperture.  If  in  space,  then 
adaptive  optics  would  not  be  needed  for  the  small  telescope.  A 
beamsplitter  in  the  small  telescope  allows  for  the  detection  of 
intensity  simultaneously  in  two  planes:  the  usual  focal  plane,  where 
there  exists  a  diffraction-limited  image  of  the  object,  and  a 
demagnified  image  of  the  aperture  plane.  The  diffraction-limited  image 
of  the  object  has  low  resolution  since  it  comes  from  a  small  aperture. 
It  is  assumed  that  the  intensity  measurements  are  made  over  a  short 
e- ough  time  that  the  object  and  the  receiver  are  essentially  fixed  in 
space  relative  to  one  another. 
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In  summary,  the  optical  receiver  makes  the  following  intensity 

measurements.  Letting  (a  binary  function)  denote  the  entire  large 

aperture  (including  the  small  aperture)  and  A  denote  the  small 

2  ^ 

aperture,  we  have  I F (u) !  [A.  (u)  -  A  (u)]  from  the  light-bucket 

2  L  o 

detectors,  IF(u)p  A  (u)  from  the  reimaged  aperture  of  the  small 

s2 

telescope,  and  I g (x) I  ,  the  intensity  of  the  low-resolution  image, 
where  g(x)  =  a$(x)  *  f(x),  ag(x)  is  the  Fourier  transform  of  A  (u)  and 
*  denotes  convolution. 

4.4  DATA  PROCESSING  STEPS 

Figure  4-2  shows  a  block  diagram  depicting  how  these  three 
intensity  measurements  (or  their  square  roots,  the  moduli  or 
magnitudes)  are  used  to  retrieve  the  phase  over  the  large  aperture  and 
reconstruct  a  fine-resolution  image.  A  support  constraint  for  the 
object  is  computed  (in  either  of  two  ways),  and  the  phase  over  the 
small  aperture  is  retrieved.  Then  the  fine-resolution  image  is 
reconstructed  using  all  the  available  information.  In  what  follows 
each  of  these  steps  is  described  in  some  detail. 

4.4.1  Support  Estimation 

A  support  constraint  for  the  object  can  be  obtained  in  one  of  two 
ways:  from  the  low-resolution  image  or  using  a  triple-intersection  of 
the  autocorrelation  support  [4.3,4.18]. 

An  estimate  of  the  support  of  the  object  can  be  obtained  from  the 

2 

low-resolutirn  image,  I  g  ( x )  i  ,  by  thresholding  it  at  an  appropriate 

o 

level  (i.e.,  the  support  function  is  set  equal  to  unity  where  l  g  ( x )  I 
exceeds  the  threshold  and  xero  elsewhere).  If  the  threshold  level  is 
set  too  high,  then  the  support  of  the  object  is  underestimated.  If  the 
tnreshold  :e->  1  is  set  too  low,  then  the  support  is  overestimated 
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Figure  4-2.  Data  processing  steps  to  reconstruct  a  fine-resolution 

imag^  (retrieve  the  phase  in  the  aperture  plane)  from  the 
intensity  measurements. 
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because  of  noise  and  sidelobes  due  to  diffraction.  The  sidelobes  can 
be  minimized  by  placing  an  apodization  (weighting)  across  the  small 
aperture,  but  at  the  expense  of  optical  efficiency  and  resolution  of 
the  'ow-resol ution  image. 

The  method  we  find  useful  for  selecting  the  threshold  level  is  as 
follows.  Several  candidate  thresholded  low-resolution  images  are 
computed  using  different  threshold  values.  When  the  threshold  value  is 
too  high,  small  changes  in  the  threshold  value  tend  to  cause  small 
changes  in  the  area  of  the  thresholded  image.  When  the  threshold  falls 
below  the  value  needed  to  pick  up  the  noise  and/or  sidelobes,  then  the 
area  of  the  thresholded  image  grows  very  rapidly,  spreading  over  the 
entire  image  space,  and  the  thresholded  image  breaks  up  very  rapidly. 
This  is  illustrated  by  the  example  shown  in  Figure  4-3  for  the  case  of 
a  telescope  with  a  small  circular  aperture  (a  diameter-16  circle 
imbedded  in  a  128  x  128  array).  Thus  a  good  choice  of  threshold  value 
is  one  just  larger  than  the  values  for  which  the  area  of  the 
thresholded  image  grows  rapidly.  Figure  4-4  shows  the  corresponding 
results  for  the  case  of  a  weighting  function  across  the  small  aperture. 
The  weighting  function  was  chosen  to  be  the  autocorrelation  of  a  circle 
of  half  the  diameter  of  the  small  aperture;  that  is,  the  weighting 
function  falls  to  zero  at  the  edges  of  the  small  aperture.  With  the 
aperture  weighting  included,  the  diffraction  sidelobes  are  greatly 
reduced  and  the  area  of  the  thresholded  image  is  much  less  sensitive  to 
changes  in  the  threshold  value,  making  the  aperture  weighting 
worthwh.ie  despite  the  loss  of  resolution  it  causes. 

For  the  case  of  diffusely  scattering  objects,  the  Fourier  intensity 
is  a  speckle  pattern  and  the  image  (the  low-resolution  image  as  well  as 
"die  fine-resolution  image),  is  speckled  [4.19],  as  can  be  seen  in 
figures  4-3  and  4-4.  ‘lulls  in  the  thresholded  image  due  to  speckle 

o 

nulls  in  I g ( x , y ) 1 ^  can  be  eliminated  by  convolving  the  thresholded 


6b 


?>ERIM 


THRESHOLDED 


RESOLUTION 


APERTURE  WEIGHTING 


RESOLUTION 


£ 


igure  4-3. 


Thresholding  the  low-resolution  intensity  image  to 
estimate  a  support  constraint,  with  no  weighting  on  the 
small  aperture.  (A)  Diffraction-limited  low-resolution 
image  (overexposed  in  order  to  show  the  sidelobes  that 
extend  beyond  the  support  of  the  object);  (B)-(D) 
thresholded  images,  with  threshold  values  equal  to  (B) 
0.073,  (C)  0.157,  and  (D)  0.392  of  the  maximum  value  of 
the  image. 


4-4.  Thresholding  the  low-resolution  intensity  image  to 

estimate  a  support  constraint,  with  a  weighting  on  the 
small  aperture.  (A)  Diffraction-limited  low-resolution 
image  (not  overexposed);  (B)-(D)  thresholded  images,  with 
threshold  values  equal  to  (B)  0.078,  (C)  0.157,  and  (D) 
0.392  of  the  maximum  value  of  the  image. 
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image  with  a  small  circle,  then  rethresholding,  as  illustrated  in 
Figure  4-5.  The  estimate  of  the  support  of  the  object,  shown  in  Figure 
4-5 (C) ,  is  used  as  a  support  constraint  in  the  final  step  of  fine- 
resolution  image  reconstruction  by  the  iterative  Fourier  transform 
algorithm.  Since  this  support  constraint  is  only  approximate  and  may 
be  too  small,  it  is  often  useful  to  enlarge  the  support  constraint  to 
ensure  that  the  object  fits  within  it.  We  typically  enlarge  the 
support  constraint  by  adding  pixels  to  the  edges  of  the  initial  support 
constraint,  as  shown  by  the  example  in  Figure  4-5(D). 

A  second  method  of  generating  a  support  constraint,  which  uses  the 
Fourier  intensity  over  the  entire  aperture,  is  the  method  of  triple 
intersection  of  the  autocorrelation  support  [4.3].  First  the  Fourier 
intensities  over  the  small  aperture  and  from  the  light-bucket  detectors 
are  combined  to  arrive  at  the  intensity  over  the  entire  aperture.  This 
intensity  is  inverse  Fourier  transformed  to  obtain  a  fine-resolution 
autocorrelation  of  the  object.  The  autocorrelation  function  is 
thresholded  and  the  nulls  due  to  speckles  are  eliminated  in  a  way 
similar  to  that  shown  in  Figure  4-5.  Noise  and  sidelobes  outside  the 
true  autocorrelation  support  are  eliminated  to  the  extent  possible  by  a 
similar  method  operating  on  the  complement  of  the  support.  The  result 
is  an  estimate  of  the  support  of  the  autocorrelation.  Then  three 
appropriate  translates  of  the  autocorrelation  support  are  intersected 
to  arrive  at  an  upper  bound  on  the  support  of  tne  object  [4.3].  In 
this  case  the  support  constraint  is  not  an  estimate  of  the  support  of 
the  object,  but  is  an  upper  bound  that  contains  all  possible  object 
supports  consistent  with  the  support  of  the  autocorrelation. 

The  support  constraint  computed  from  the  autocorrelation  function 
is  from  fi ner- resol ut ion  data  and  therefore  may  be  more  accurate,  but 
it  may  also  be  too  large  since  reconstruction  of  the  support  of  an 
object  from  the  support  of  its  autocorrel ation  function  is  ambiguous 
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Figure  4-5.  Removal  of  nulls  due  to  speckles  in  the  image.  (A) 

Thresholded  image  from  Figure  4-4  (C) ,  (B)  convolution  of 
(A)  with  a  circle  of  diameter  7  pixels  (about  half  the 
diameter  of  a  speckle),  (C)  thresholding  of  (B)  at  0.58  of 
its  peak,  (D)  enlarged  version  of  (C)  that  may  be  used  to 
ensure  that  the  object  fits  within  it. 
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for  wide  classes  of  objects  [4.20,4.3].  The  thresholded  low-resolution 
image  avoids  these  ambiguity  problems,  but,  being  lower  resolution,  may 
not  be  as  accurate.  Further  study  is  required  to  determine  which  of 
the  two  methods  is  better  and  to  devise  a  way  to  combine  the  best 
features  of  each  into  a  composite  support  estimate. 


4.4.2  Smal 1 -Aperture  Phase  Retrieval  Using  Gerchberg-Saxton 


The  phase,  'i(u),  of  the  optical  field,  F(u)  As(u),  in  the  plane  of 
the  small  aperture  is  determined  from  the  intensities  in  the  focal 
plane  and  the  image  of  the  aperture  plane  using  a  variation  of  the 
Gerchbe  q-Saxton  algorithm  that  is  accelerated.  Figure  4-6  shows  a 
block  diagram  for  the  Gerchberg-Saxton  algorithms.  Here  we  will  refer 
to  the  original  Gerchberg-Saxton  algorithm  [4.11,4.12]  as  GS  and  the 
accelerated  versions  as  GS1  and  GS2  [4.13-4.15],  the  latter  having  the 
image-domain  operation  [4.14,  Eqs.  (9)  and  (10)] 


Wx> 


gk<x> 


*  P 


2  ig(x)  I 


gk(x) 

Tiqrmr  -  gk(x>  - lg(x)l 


gdx> ' 
igk(x) i 


(4-1) 


where  p  is  a  constant,  I g ( x) I  is  the  modulus  of  the  low-resolution 
image,  g^(x)  is  the  input  image  to  the  kth  iteration  and  g ' ^ (x)  is  the 
output  image  from  the  kth  iteration.  The  rates  of  convergence  for 
these  three  algorithms  were  compared  and  p  was  optimized.  The 
differences  in  the  convergence  rates  are  affected  not  only  by  the 
choice  of  p  but  also  by  the  choice  of  the  random  phase  used  as  the 
initial  estimate.  It  was  found  that  GS2  generally  converges  much 
faster  than  GS1  which  in  turn  converges  significantly  faster  than  GS. 
A  better  method  than  using  any  single  algorithm  is  to  combine  GS2  and 
GS:  perform  several  iterations  with  GS2,  which  initially  converges 
quickly,  then  finish  off  with  several  iterations  of  GS,  which  is  more 
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Figure  4-6.  Block  diagram  of  the  Gerchberg-Saxton  algorithms.  The 
object  domain  constraint  is  the  square  root  of  the 
measured  intensity  of  the  low-resolution  image  and  the 
Fourier  constraint  is  the  square  root  of  the  measured 
intensity  over  the  small  aperture. 
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stable  and  converges  to  a  smaller  error.  For  these  algorithms  the 
object  domain  error  metric  (ODEM) ,  a  normalized  root-mean-squared  (RMS) 
error,  is  given  by 


ODEM2 


z 


2 


ig£(x)i  -  ! g (x) i 


Z  ig(x)i2 

X 


(4-2) 


which  is  a  measure  of  how  close  the  output  image  modulus  agrees  with 
the  modulus  of  the  measured  low-resolution  image  and  is  the  criterion 
by  which  we  judge  whether  the  algorithm  has  converged.  (A  similar 
error  metric  in  the  Fourier  domain  can  also  be  used.)  Note  that  in 
orcer  for  Eq.  (4-2)  to  be  meaningful  it  is  necessary  to  normalize 
I g (x , y ) I  so  that  it  has  the  same  energy  (sum  of  squares)  as  the  Fourier 
modulus  data.  The  quality  measure  we  use  to  evaluate  the 
reconstruction  results  is  the  absolute  error  of  the  complex-valued 
reconstructed  image,  also  a  normalized  RMS  error,  which  is  given  by 


Z  I*  g' (*  -  xQ)  -  g (x) 1 2 

ABSERR2  =  — - 5 - 

Z  Ig(x)i 


(4-3) 


where  a  is  the  complex  factor,  and  xq  is  the  shift,  that  minimizes 
ABSERR.  It  can  be  shown  that  xQ  is  given  by  the  location  of  the 
maximum  magnitude  of  rg'g(x)>  the  cross-correlation  of  g1  with  g;  and  a 


This  absolute  error  can  only  be  computed  in 


=  r?.g(x0)/Elg(x)l‘ 

digital  simulation  experiments  for  which  the  true  image  is  known. 
Although  it  measures  the  error  in  the  complex  numbers,  which  includes 
both  magnitude  and  phase  errors,  ABSERR  correlates  well  with  the 
standard  deviation,  o of  the  error  of  the  phase  retrieved  over  the 
small  aperture.  As  shown  in  the  Appendix,  the  expected  relationship, 
ignoring  errors  in  the  modulus,  is 
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ABSERR2  ~  1  -  exp[-cr^2]  .  (4-4) 

Figure  4->  shows  examples  of  algorithm  convergence.  Twenty 
iterations  of  either  GS2  or  GS  were  followed  by  20  iterations  of  GS. 
The  optimum  value  of  p  was  found  to  be  about  1.5  to  2.  The  algorithm 
is  not  very  sensitive  to  small  changes  in  the  value  of  p.  Retrieval  of 
the  phase  over  the  small  aperture  was  found  to  be  relatively  fast  (only 
about  30  iterations  are  required). 

To  test  the  sensitivity  of  the  combined  algorithm  to  noise,  low 
light  levels  (quantum-limited  measurements)  were  simulated  by 
subjecting  the  intensity  measurements  in  both  planes  to  a  Poisson  noise 
process.  We  chose  to  simulate  the  same  number  of  photons  in  each  of 
the  two  planes.  After  scaling  the  intensity  data  to  have  a  given 
expected  total  number  of  detected  photons,  each  pixel  was  replaced  with 
a  sample  drawn  from  a  Poisson  distribution  with  mean  and  variance  equal 
to  the  pixel  value. 

For  these  experiments,  the  object  is  approximately  of  size  40  x  60 

pixels  imbedded  in  a  128  x  128  array.  Therefore  the  intensity  of  the 

Fourier  transform  of  the  object  (computed  using  an  FFT)  is  a  speckle 

pattern  with  about  3x2  samples  per  speckle.  The  Fourier  data  was  set 

to  zero  outside  a  circle  of  diameter  16  pixels  to  simulate  the  effect 

of  the  small  aperture  (without  weighting).  Therefore  there  should  be 
2 

about  t8  / 6  ^  33  speckles  in  the  small  aperture. 

Figure  4-8  shows  the  convergence  of  the  algorithm  for  a  variety  of 
noise  levels.  Figure  4-9  shows  ABSERR,  the  quality  of  the  output 
image,  as  a  function  of  iteration  number  for  a  variety  of  noise  levels. 
Figure  4-10  shows  ABSERR  for  the  reconstructed  image  as  a  function  of 

the  total  number  of  detected  photons.  Very  good  results  are  obtained 

4  4 

for  10  or  more  photons,  which  corresponds  to  10  / 3 3  ~  300  photons  per 
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Figure  4-9.  RMS  error  (ABSERR)  of  the  complex-valued  reconstructed 

low-resolution  image  as  a  function  of  iteration  number  for 
a  variety  of  light  levels. 
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speckle.  From  these  results  we  see  that  the  Gerchberg-Saxton 
algorithms  converge  rapidly  and  are  reasonably  robust  in  the  presence 
of  noise. 

Comparing  Figures  4-8  and  4-9,  one  sees  that  ODEM  >  ABSERR,  which 
is  contrary  to  what  one  would  ordinarily  expect.  This  can  be  explained 
as  follows.  In  our  simulation  the  Fourier-  (aperture-)  domain  data  was 
only  slightly  oversampled  whereas  the  image-domain  data  was  highly 
oversampled,  that  is,  there  are  many  pixels  per  speckle.  Therefore, 
even  though  both  planes  have  the  same  average  number  of  photons  per 
speckle  (the  number  of  speckles  being  the  same  in  both  domains),  the 
’mage  domain  received  far  fewer  photons  per  pixel.  This  causes  the 
image-domain  data  to  have  a  much  greater  mean-squared  error  than  the 
Fourier-domain  data.  Therefore,  the  output  image,  g'(x),  which  is 
consistent  with  the  Fourier  modulus  data,  has  an  error  appropriate  to 
the  lower  error  of  the  Fourier  modulus  data,  and  it  is  closer  to  the 
true  image  (as  measured  by  ABSERR)  than  it  is  to  the  noisy  image-plane 
modulus  data  (as  measured  by  ODEM). 

Since  the  image  domain  was  highly  oversampled,  we  performed  a 
simple  noise  filtering.  The  noisy  image  intensity  was  Fourier 
transformed,  the  Fourier  transform  was  set  to  zero  outside  a  circle 
diameter  32  pixels  (since  the  ideal  complex-valued  image  has  a  Fourier 
transform  that  is  zero  outside  a  circle  of  diameter  16  pixels),  and  the 
result  was  inverse  Fourier  transformed  to  yield  a  smoothed  image  with 
•"educed  noise.  Before  taking  the  square  root  to  compute  the  image 
modulus,  small  negative  numbers  introduced  by  the  filtering  process 
were  set  to  zero.  Figures  4-11  and  4-12  show  the  convergence  and  image 
quality  for  a  set  of  experiments  similar  to  the  ones  described  above, 
but  using  the  filtered  images.  As  expected,  ODEM  was  lower  for  the 
case  of  noise  filtering  than  without  it.  However,  ABSERR,  which  is 

ultimately  of  greater  importance,  was  slightly  better  without 
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:er'ing;  consequently  it  is  better  net  tc  ’  mvi  s 

• •  wee.  This  filtering  issue  require'  - 

tote  from  Figures  4-8  and  4-11  that  the  initial  convergence  rate  of 
■  i.gorithm  is  not  monotonical ly  related  to  the  amount  of  noise 
merit.  This  results  from  the  fact  that  some  of  the  ranee m  phases 
: e  J  as  the  initial  estimate  are  closer  to  me  true  m 1  mm  .nan 
- h  e  rs  . 


Fine-Resolution  Image  Reconstruction: 

Wi  th  all  the  data  in  hand  --  including  the  Fourier  mmrumv  over 
m  entire  aperture,  the  phase  over  the  sma :  1  ,:pert..re;  and  mm  i  up  port 
,'mtraint  --  we  perform  image  reconstruction  mo..  ■  -  o;vp 

-•  mr  transform  algorithm,  which  seeks  a  solution  mr-mm---'  ■  .-.'ll 

e  data  and  constraints  [4.1,  4. 7-4.9,  4 . 1 3  -  ► .  1  :  ’  .  o  ■  -  m-mm  of 
-  algorithm,  which  is  a  generalization  : f  the  urn  ■  ■  t.::i 

-gorithm,  is  shown  in  Figure  4-13. 

When  using  any  phase  information  in  me  Four1  ^  the 

m  m  ion  of  the  support  constraint  in  the  •  •  -  a : ;  m.  chosen 

.  ..  :m  consi stent  with  the  given  Fourier  pnase.  { F o r  ..  .  ntional 

;  o  retrieval  with  no  a  priori  phase  information,  it  dees  not 

‘  m)  One  way  to  ensure  this  is  to  crcss-co- ’  a m  :  ue  saocort 

traint  with  the  low-resolution  image  and  use  the  me  mm  n  m  the 
value  of  the  cross-correlation  to  determine  the  optimal  position 
of  the  support  constraint. 

For  the  case  of  a  di  f f icul  t-to-reco.istruct  comp!  ex- val  ued  object, 
initial  attempts  to  use  the  smal 1 -aperture  phase  with  the  iterative 
Fourier  transform  algorithm  were  unsuccessful,  whether  the  phase  was 
just  used  in  the  initial  estimate  or  reinforced  during  the  iterations. 
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igure  4-13.  Block  diagram  of  the  iterative  Fourier  transform 

algorithm.  The  object-domain  constraint  is  a  support 
constraint  derived  from  the  measured  data,  and  the 
Fouri er-domai n  constraints  are  the  square  root  of  the 
measured  intensity  over  the  entire  large  aperture  and  the 
phase  retrieved  by  the  Gerchberg-Saxton  algorithm  over 
the  small  aperture. 
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A  possible  reason  for  this  failure  is  the  fact  that  the  area  of  the 
small  aperture  is  a  small  fraction  of  the  area  of  the  entire  large 
aperture,  and  so  the  incorrect  phase  over  the  rest  of  the  aperture 
overwhelms  the  influence  due  to  the  correct  phase  over  the  small 
aperture. 

The  foPowing  modifications  to  the  algorithm  were  found  to  be 
necessary  for  a  reliable  reconstruction.  First,  in  order  to  reduce 
impulse-response  sidelobes  it  is  advantageous  to  use  a  weighting 
function  on  the  Fourier  modulus.  This  is  important  since  large 
sidelobes  extending  beyond  the  edges  of  the  object  support  will  violate 
the  support  constraint  and  hinder  convergence.  The  Fourier  modulus 
weighting  function  was  chosen  to  be  the  autocorrelation  of  a  circle. 
Initially  the  diameter  of  the  circle  was  chosen  to  be  such  that  the 
Fourier  modulus  weighting  went  to  zero  over  an  area  just  slightly 
larger  than  the  area  over  which  the  small -aperture  phase  was  known. 
Then  a  cycle  of  30  hybrid  input-output  iterations  followed  by  10  error- 
reduction  iterations  [4.1]  was  performed,  reinforcing  the  phase  over 
the  small  aperture  at  each  iteration.  (In  practice,  reasonably  good 
results  can  also  be  obtained  if  the  phase  over  the  small  aperture  is 
used  for  the  first  iteration  only,  without  being  reinforced  during 
later  iterations;  but  better  results  are  obtained  by  continual 
reinforcement  of  the  known  phase.)  This  was  sufficient  to  converge  to 
a  solution  for  the  phase  over  the  nonzero  area  of  the  weighted  Fourier 
modulus,  since  the  phase  was  already  known  over  most  of  that  area  to 
begin  with.  Then  the  Fourier  modulus  was  re-weighted  with  a  weighting 
function  of  slightly  larger  area,  and  another  cycle  of  iterations  was 
performed.  This  process  was  continued  until  the  weighting  function 
encompassed  the  entire  area  of  the  measured  Fourier  modulus  data.  Thus 
the  phase  retrieval  proceeded  by  a  bootstrap  approach,  with 
successively  larger  areas  of  phase  retrieved,  and  successively  finer- 
resolution  images  reconstructed,  during  each  cycle  of  iterations.  When 


83 


Srim 


we  compute  ABSERR  by  Eq.  (4-3),  we  use  for  g(x)  the  diffraction-limited 
image  for  the  same  weighting  of  the  Fourier  transform  as  is  being  used 
for  the  Fourier  modulus  weighting  for  that  cycle  of  iterations.  We 
have  recently  learned  that  others  have  also  found  an  expanding  weighted 
modulus  approach  to  be  important  for  reconstructing  complex-valued 
images  [4.21]. 

Ordinarily  when  reconstructions  are  performed  using  a  poorly-known 
support  constraint  (as  is  the  case  here),  we  have  found  it  best  to 
start  with  a  smaller  support  constraint  for  early  iterations  and  expand 
the  support  constraint  for  later  iterations.  However,  when  using  the 
expanding  weighted  modulus  algorithm,  the  images  that  are  reconstructed 
during  the  early  iterations  are  larger  than  the  images  reconstructed 
during  the  later  iterations,  since  for  the  early  iterations  the  point- 
spread  function  is  much  larger  due  to  the  use  of  a  narrow  weighting 
function  in  the  Fourier  domain.  By  experimentation  with  support 
constraints  that  were  expanded  or  shrunk  as  the  iterations  progressed, 
we  found  that  a  good  strategy  was  to  use  a  support  constraint 
appropriate  for  the  low-resolution  image,  and  keep  it  fixed  during  all 
the  iterations.  However,  an  alternate  strategy  may  be  necessary 
depending  on  the  ratio  of  the  diameters  of  the  small  and  large 
apertures  or  on  how  the  support  constraint  is  formed. 

For  fine-resolution  image  reconstruction  from  the  Courier  modulus, 
for  which  the  only  image-domain  constraint  is  a  support  constraint,  the 
object-domain  error  metric  is  given  by,  instead  of  Eq.  (4-2), 


H  1 (x) I  2 

ODEM2  =  ^ x 

H  ! (x) ! 


i.e.,  the  energy  outside  the  support  constraint  S. 


(4-5) 
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4.5  IMAGE  RECONSTRUCTION  EXAMPLE 

Figure  4-14  shows  an  example  of  image  reconstruction  using  the 
approach  described  above.  In  the  upper  left  is  the  Fourier  modulus 
data  (noise  free)  over  the  entire  aperture,  a  circle  of  diameter  64 
pixels  (imbedded  in  a  128  x  128  array)  with  the  aperture  of  the  small 
telescope  indicated  by  a  dark  circle.  In  the  upper  right  is  the  low- 
resolution  image  obtained  through  the  small  aperture  of  diameter  16 
pixels,  weighted  by  the  autocorrelation  of  a  diameter-8  circle.  The 
support  constraint,  shown  in  the  lower  left,  was  obtained  by 
thresholding  the  1 ow-resol uti on  image  as  described  above,  and  the 
smal 1 -aperture  phase  was  estimated  using  the  accelerated  Gerchberg- 
Saxton  algorithm.  All  that  information  --  the  Fourier  modulus  over  the 
large  aperture,  the  support  constraint,  and  the  Fourier  phase  over  the 
small  aperture  --  was  combined  to  retrieve  the  phase  over  the  large 
aperture  by  the  iterative  Fourier  transform  algorithm  using  the 
expanding  weighted  modulus  approach.  After  25  cycles  of  iterations 
during  which  the  weighting  was  expanded,  plus  an  additional  6  cycles  at 
the  end,  for  a  total  of  over  one  thousand  iterations,  the  image  shown 
in  the  lower  center  was  obtained.  It  is  very  close  to  the  true  fine- 
resolution  diffraction-limited  image,  shown  in  the  lower  right.  Figure 
4-15  shows  intermediate  results  with  different  weightings  on  the 
Fourier  modulus.  The  phase  is  retrieved  wll  for  each  weighting  of  the 
Fourier  modulus  before  the  weighting  lunction  is  expanded,  and  so  at 
each  step  a  diffraction-limited  image  (for  the  resolution  given  by  the 
weighting  function)  is  reconstructed .  Figure  4-16  shows  the  object- 
domain  error  metric  (ODEM)  and  the  absolute  error  (ABSERR)  as  a 
function  of  iteration  number.  Also  indicated  is  the  diameter  of  the 
Four! er-wei ghti ng  function  as  the  iterations  progress.  It  was  found 
that  if  substantially  fewer  iterations  per  cycle  were  used  or  if  larger 
jumps  in  the  size  of  the  weighting  function  were  used,  then  the 
convergence  of  the  algorithm  was  much  less  reliable.  Unlike  previously 
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Figure  4-14.  Image  reconstruction  example.  (Left  to  right.)  (A) 

Fourier  modulus  data  over  a  large  circular  aperture  -- 
the  black  circle  shows  the  area  of  the  small  apert  .  e; 

(B)  low-resolution  image  from  the  small  aperture; 
object  support  constraint  derived  from  (B);  (D)  image 
reconstructed  by  the  Gerchberg-Saxton  algorithm  followed 
by  the  iterative  Fourier  transform  algorithm  using  (A). 
(B)  and  (C);  (E)  ideal  image  for  comparison. 
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IfWGE  RECONSTRUCTION  USING  PHASE  OUER  SHALL  APERTURE 


<B>  21  REF  <0>  31  REF  <F)  43  REF  <H>  63  REF 


<C>  ERR  -  .063  CE>  ERR  -  .987  <G>  ERR  -  .069  <1>  ERR  »  .122 


Figure  4-15.  Intermediate  reconstruction  results  w'th  different 
weightings  on  the  Fourier  modulus.  Bottom  row: 
reconstructed  images;  top  row:  ideal  images  with  the 
same  Fourier  weighting.  Diameter  of  weighting  function 
in  pixels:  (A)  21,  (B)  31.  (C)  43,  {O'  53. 
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Hgijro  4-15.  Convergence  of  the  iterative  Fourier  transform  algorithm: 

ODEM  (□)  and  ABSERR  (A)  (solid  lines  --  scale  at  left) 
and  the  diameter  of  the  Fouri er-modul us  weighting 
function  (dashed  line  --  scale  at  right)  as  a  function  of 
iteration  number. 
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published  results,  in  which  ODEM  (or  ABSERR)  starts  out  large  and 
decreases  with  iteration  number,  here  it  starts  low  and  stays  low  since 
at  any  given  point  we  are  trying  to  retrieve  only  an  additional  thin 
annulus  of  phase.  A  sawtooth  behavior  is  seen  since  the  error  jumps  up 
each  time  the  weighting  on  the  Fourier  modulus  is  enlarged. 

4.6  CONCLUSION 

A  complex-valued  image  of  an  object  with  convex  support  and  without 
bright  glints,  whose  support  is  not  known  a  priori,  is  ordinarily  very 
difficult  to  reconstruct  from  its  Fourier  modulus.  We  have 
demonstrated  that  a  low-resolution  intensity  image  of  the  object,  taken 
through  a  smal 1 -aperture  telescope  contiguous  with  intensity 
measurements  over  a  large  aperture,  can  be  used  to  help  reconstruct  a 
fine-resolution  image.  The  low-resolution  image  is  used  both  to 
determine  the  Fourier  phase  over  the  small  aperture  and  to  form  a 
support  constraint  for  the  object.  The  retrieval  of  the  phase  over  the 
small  aperture,  using  an  accelerated  version  of  the  Gerchberg-Saxton 
algorithm,  was  found  to  be  not  only  fast,  but  also  robust  in  the 
presence  of  noise.  The  reconstruction  of  the  fine-resolution  image  was 
also  successful,  but  was  found  to  take  a  large  number  of  iterations. 
The  determination  of  its  performance  in  the  presence  of  noise  will 
require  further  research,  which  is  presently  being  planned. 

Portions  of  this  chapter  were  presented  at  the  O.S.A.  Topical 
Meeting  on  Signal  Recovery  and  Synthesis  III.  N.  Falmouth,  MA,  14-16 


June  1989  [4.22]. 
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5 

IMAGE  RECONSTRUCTION  LABORATORY  EXPERIMENTS 

5.1  SUMMARY 

A  new  method  for  image  reconstruction  from  Fourier  intensity  and 
low-resolution  image  intensity  measurements  was  introduced  in  Chapter  4 
and  computer  simulation  experiments  to  demonstrate  this  method  were 
described.  An  important  next  step  in  the  development  of  this  method  is 
the  investigation  of  image  reconstruction  from  experimental  data 
collected  in  the  laboratory.  Therefore,  an  experiment  was  set  up, 
laboratory  data  was  collected,  and  image  reconstruction  procedures  were 
applied  with  and  without  the  use  of  the  low-resolution  image.  An  image 
was  reconstructed  using  Fourier  intensity  measurements  alone.  Further 
work  to  improve  data  collection  procedures  is  required  to  reconstruct 
an  image  using  the  additional  low-resolution  image. 

The  remainder  of  this  section  gives  further  details.  In  Section 
5.2,  it  is  shown  that  accurate  Fourier  intensity  data  can  be  collected 
in  the  Fresnel  zone  of  the  illuminated  test  object.  In  Section  5.3, 
experimental  work  to  verify  the  linearity  and  spatial  uniformity  of  the 
2-D  CCD  detector  is  discussed.  The  experimental  set-up  is  described  in 
Section  5.4.  Necessary  data  sampling  and  pre-processing  is  discussed 
in  Section  5.5  and  image  reconstruction  results  are  given  in  Section 
5.6. 

5.2  FRESNEL  ZONE  DATA  COLLECTION 

Image  reconstruction  procedures  based  on  phase  retrieval  require 
measurement  of  the  Fourier  intensity  of  the  illuminated  object.  It  is 
well-known  that  the  Fourier  intensity  can  be  measured  in  the  Fraunhofer 
region  with  respect  to  the  object  [5.1].  However,  for  imaging  of  space 
objects  from  space-based  sensors,  the  sensor  is  more  likely  to  be  in 
the  Fresnel  region  than  the  Fraunhofer  region  with  respect  to  the 


93 


object.  For  example,  for  an  object  of  diameter  d  =  1  m,  illuminated  by 
a  laser  of  wavelength  X  =  1  and  a  sensor  at  a  distance  z  with 
aperture  diameter  D  =  1  m,  the  Fraunhofer  approximation  requires  [5.1] 

z  »  Td2/4X  :  1  Mm  (5-1) 
while  the  Fresnel  approximation  requires  [5.1] 

z3  »  ir(d  +  D)4/64X  =  (100  m)3  .  (5-2) 

A  sensor  at  a  likely  distance  of  z  =  1  Mm  would  therefore  easily 
satisfy  the  Fresnel,  but  not  quite  satisfy  the  Fraunhofer, 
approximation.  A  similar  situation  exists  in  a  laboratory  experiment 
unless  a  Fourier  transform  lens  is  used  between  the  object  and  the 
detector. 


When  the  Fresnel  approximation  is  satisfied,  the  relationship 
between  the  optical  field  U^x,  y)  resulting  from  illumination  of  the 
object  and  the  field  U2(u,  v)  at  a  distance  z  is  [5.1] 


,(U,  V)  .  expp'^(z-  -2)]  /  /  V*.  y> 


2  ± 


exp 


-i2r(xu  +  yv 


Xz 


dx  dy  . 


(5-3) 


The  intensity  I(u,  v)  is 


I(u,  v)  =  IU2(u,  v)l2  =  (Xz)"2  J  f  U 1 ( x ,  y)  exp 


1t(x2  i  y2) 
Xz 


_  f-i 2t(xu  +  yv) 

*  exp[ - xT 


dx  dy 


(5-4) 
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If  phase  retrieval  is  successfully  applied  to  measurements  of  the 
Fourier  intensity  data  I(u,  v) ,  the  complex-valued  image  reconstructed 
will  be  the  di ffraction-1 imi ;ed  image  of  U^(x,  y)  exp[i*-(x  +  y  )/Xz], 
whereas  a  diffraction-limited  image  of  U^(x,  y)  is  desired.  (The 
diffraction  limit  will  be  determined  by  the  aperture  diameter  D  over 
which  Fourier  intensity  data  is  collected.) 

2  2 

However,  the  image  obtained  of  U^(x,  y)  exp[ir(x  +  y  )/Xz]  can  be 
equally  useful  to  that  of  U^(x,  y)  under  the  following  conditions.  If 
the  object  is  rough  compared  to  the  wavelength  of  the  illumination, 
then  the  quadratic  phase  is  added  to  the  random  phase  of  U,(x,  y).  The 

-*•  9 

intensities  of  the  diffraction-limited  images  of  U,(x,  v)  exp[i*-(x"  + 
2  i 

y  )/Xz]  and  U^x,  y)  are  then  different  realizations  (speckle  patterns) 

of  the  same  image  of  the  object.  If,  in  addition,  it  is  not  desired  to 

use  the  phase  information  of  U^x,  y),  then  either  image  gives  the  same 

information  and  reconstruction  of  images  from  intensity  data  measured 

in  the  Fresnel  region  of  the  object  is  as  acceptable  as  that  from  data 

measured  in  the  Fraunhofer  region. 

5.3  DETECTOR  CALIBRATION 

It  is  known  from  previous  work  at  ERIM  that  errors  in  the  Fourier 
intensity  data  lead  to  errors  in  reconstructed  images  and,  if  the  error 
is  too  large,  to  an  inability  to  reconstruct  an  image  [5.2,  5.3].  Both 
the  optical  propagation  path  from  the  object  to  the  detector  and  the 
detector  itself  can  cause  data  errors.  Propagation  path  errors  depend 
on  the  experimental  set-up  and  are  discussed  in  Section  5.4.  Possible 
errors  which  can  be  introduced  by  a  detector  include  nonlinear  response 
to  input  intensity,  spatially  nonuniform  response,  and  additive  noise. 
These  detector  errors  are  the  subject  of  this  section. 

A  Fairchild  CCD3000F  camera  using  the  Fairchild  CCD222  sensor  with 
380  by  488  pixels  had  been  used  in  earlier  successful  phase  retrieval 
demonstration  experiments  [5.4]  and  was  chosen  for  the  current 
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experiments  as  well.  The  linearity  of  the  response  was  carefully 
measured  using  a  controllable  light  source  and  calibrated  detector. 
Spatial  uniformity  of  response  was  also  measured  in  this  way  for  all 
pixels.  The  linearity  of  the  typical  CCD  pixel  was  ±2%  and  the  spatial 
uniformity  (standard  deviation  computed  over  all  pixels)  was  1.5%  when 
measured  at  a  light  level  corresponding  to  the  maximum  signal 
obtainable  when  still  operating  in  the  quasilinear  response  region. 
Software  was  developed  so  that  the  response  of  each  pixel  could  be 
independently  corrected  for  nonlinearity  and  bias.  However,  the 
correction  was  found  to  be  so  small  that  it  did  not  affect  the  results. 
The  standard  deviation  of  the  frame-to-frame  noise  of  the  Fairchild  CCD 
camera  was  measured  and  found  to  be  <  1%,  again  measured  at  the  light 
level  corresponding  to  the  maximum  quasilinear  response. 

5.4  FXPERIMENTAL  OPTICAL  SET-UP 

The  experimental  optical  set-up  is  shown  in  Fig.  5-1.  An  Argon  ion 
laser  operating  at  a  wavelength  X  =  0.5145  pm  was  used  to  illuminate  a 

transmissive  object  consisting  of  a  binary  transmittance  mask  placed  in 
contact  with  a  ground  glass.  The  maximum  width  of  the  object  was  about 
10  mm.  The  transmitted  light  propagated  1.54  m  to  a  lens  of  focal 
length  f^  =  1.54  m.  The  lens  was  about  40  mm  thick  in  the  center.  As 
explained  later  in  this  section,  this  lens  is  used  to  cancel  the 
quadratic  phase  due  to  propagation  over  a  distance  f^.  Two  apertures 
A^  and  A2  were  located  behind  the  last  surface  of  the  lens  ol  distances 
of  8  and  22  mm,  respectively.  For  collection  of  Fourier  intensity 
data,  the  CCD  detector  was  placed  in  plane  P^  about  20  mm  behind  the 
second  aperture.  For  collection  of  image  data,  the  optical  field  at 
plane  P^  was  Fourier  transformed  by  a  lens  to  a  plane  P^  to  which 
the  CCD  detector  was  moved.  The  imaging  lens  had  a  focal  length 
=  189  mm  and  was  positioned  so  that  the  planes  P^  and  P^  corresponded 
to  its  front  and  back  focal  planes,  respectively. 
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Figure  5-1.  Optical  Setup  for  Phase  Retrieval  Laboratory  Experiments. 
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The  data  collected  consisted  of  Fourier  intensity  data  over  both  a 
]arge  and  a  small  aperture  and  image  data  corresponding  to  both  the 
large  and  small  apertures.  The  apertures  were  square  and  approximately 
5  mm  and  1.2b  mm  in  width,  respectively.  The  small  aperture  (A,,) 
Fourier  intensity  data  and  the  low-resolution  image  data  (obtained 
using  the  small  aperture)  were  input  to  the  Gerchberg-Saxton  algorithm 
described  in  Chapter  4.  For  high-resolution  image  reconstruction,  the 
large  aperture  (A^)  Fourier  intensity  data  can  be  input  to  an  iterative 
Four’er  transform  algorithm  using  in  addition  either  (1)  an  upper  bound 
on  object  support  estimated  by  the  triple  intersection  method  [5.5]  or 
(2)  an  object  support  estimated  by  thresholding  the  low-resoiution 
image.  In  the  second  case,  the  Fourier  phase  over  the  small  aperture, 
estimated  by  the  Gerchberg-Saxton  algorithm,  can  also  be  used.  The 
first  method  was  investigated  via  computer  simulations  in  earlier  work 
at  ERIM  [b.6].  The  second  method  is  the  new  image  reconstruction 
method  described  in  Chapter  4. 

Because  of  its  relative  simplicity,  the  optical  system  was  found  to 
be  sufficiently  stable  that  many  frames  of  Fourier  or  image  data  could 
be  averaged  to  reduce  additive  detectGr  noise.  The  optical  path  from 
the  object  to  the  detector  was  enclosed  in  black  velvet  to  reduce  stray 
1 i ght  to  a  jery  1 ow  1 evel  . 

The  object  diameter  d  was  about  10  mm,  the  distance  z  from  the 
object  to  the  detector  in  plane  was  about  1.5  m,  and  the  maximum 
detector  width  D  was  about  5  mm  correspond!' ng  to  the  use  of  the  central 
256  by  256  pixels  of  the  detector  array.  Substituting  these  values 
into  Fq.  (5-2)  for  the  Fresnel  approximation  gives 


z  J  »  (0.17  ml' 


(5- 5a) 
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while  the  Fraunhofer  approximation  [Eq.  (5-1)]  would  require 

z  »  150  m  .  ( 5 - 5b ) 

The  detector  is  therefore  in  the  Fresnel  region  of  the  object  and 
Fourier  intensity  data  can  be  collected  as  described  in  Section  5.2. 

It  should  be  noted  that  the  small  and  "large  apertures  are  placed  in 
front  of  the  CCD  detector  by  about  20  and  34  mm,  respectively. 
However,  the  plane  P^  where  the  Fourier  data  is  collected  and  the  plane 
P-,  where  the  image  data  is  collected  are  related  by  a  Fourier 
transformation  because  they  arc  the  front  and  back  focal  planes  of  the 
lens  . 

The  lens  is  required  to  cancel  the  quadratic  phase  [the 
exponential  in  front  of  the  integral  in  Eq.  (5-3)]  due  to  propagation 
from  the  object  to  plane  P^  so  that  the  Fourier  transformation  by  the 
lens  l_2  of  the  optical  field  in  plane  P^  gives  a  focused  image  in  the 
plane  P^.  The  image  formed  at  plane  P^  will  include  the  quadratic 
phase  described  in  Section  5.2. 

Ideally,  the  lens  and  the  two  apertures  would  be  located  in  the 
plane  P^.  The  lens  wouiu  then  not  change  the  Fourier  intensity 
measurement  since  the  lens  and  the  detector  would  be  in  the  same  plane. 
Physically,  this  is  impractical.  The  phase  error  in  plane  Pj  caused  by 
the  approximately  42  mm  distance  between  the  lens  and  the  plane  Pj 
can  be  determined  as  follows.  The  optical  field  immediately  after  the 
lens  Lj  is  given  by 


U(u,  v)  = 


exp(i2r  z/X) 


iXz 
•  exp 


f  {  Uj(x,  y) 


exp 


\r(x2  +  y2' 1 


Xz 


-i2r(xu  +  yv) 
Xz 


dx  dy 


(5-6) 
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where,  as  shown  in  Fig.  5-1,  z  -  The  angular  spectrum  A  (f  ,  f  ) 

of  this  field  is  [5.1] 

Ao(fu'  fv>  =  VXzfu-  Xzfv}  exp[iT\z(fu2  +  f/) 

where 

fu  =  x/Xz  and  fy  =  y/Xz  .  (5-7) 

The  angular  spectrum  A(f  ,  fv)  after  propagation  an  additional  distance 
Az  is 


Mfu.  V  ■  W  fv)  «p{ 


i2*  Az 


>  -  X2  fu2  -  X?  f/ 


1/2 


}  .  (5-8) 


The  field  after  this  propagation  is  therefore 


2  ,  ..2s 


''(u.  v)  .  gxR(i|,.z,A),  f  f  Ui(Xi  y)  e>p[Ut»  ..J  1  X 


•  exp 


i2r  Az 


2  2 
l  _  x  -LI 


1/2, 


exp 


-12t(xu  +  yv) 


Xz 


dx  dy. 
(5-9) 


The  field  at  the  plane  differs  from  that  after  the  lens  Lj  by  a 
phase  factor 


exp 


i2*  Az 


2  2 

1  .  x  _+  ,y 


il/2. 


(5-10) 


within  the  integral  in  Eq.  (5-9).  The  only  effect  of  this  phase  factor 
is  once  again  to  change  the  realization  of  the  image  of  the  object. 
[See  the  discussion  after  Eq.  (5-4).]  Note  that  aperture  effects  have 
not  been  included  in  this  initial  analysis. 
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Experimentally,  the  2-D  speckle  intensities  at  a  plane  behind  the 
lens  and  at  another  plane  50  mm  behind  the  first  were  detected  and 
found  to  be  very  nearly  equal.  This  indicates  that  the  42  mm 
propagation  distance  between  the  lens  and  the  plane  does  not 

adversely  affect  the  data  collected. 

5.5  DETECTOR  SAMPLING  AND  DATA  PREPROCESSING 

For  a  rectangular  diffuse  object  of  width  d  and  d  in  orthogonal 

x  y 

directions,  the  expected  value  of  the  normalized  autocorrelation  of  the 
intensity  of  the  speckle  pattern  at  distance  z  is  [5.7]: 


R(u,  v)  =  1  +  sine 


fd  u] 

fd  v] 

x 

IXz  J 

si  nc 

JL. 
Ixz  J 

(5-11) 


where  X  is  the  wavelength  of  illumination  and  u  and  v  are  orthogonal 
coordinates.  The  first  zeros  of  the  second  term  in  the  autocorrelation 
function  are  located  at  u  =  Xz/dx  and  v  =  Xz/d^.  The  Nyquist  sampling 
rates  are  half  these  distances,  so  the  detector  separations  Au,  Av  must 
sati sfy 


Au  <  ,  Av  <  .  (5-12) 
x  y 

For  the  objects  used  in  these  experiments,  the  maximum  width  was  about 
10  mm,  so  the  Nyquist  sample  spacing  was  40  fim  in  each  direction.  For 
the  Fairchild  CCD  detector,  the  center-to-center  detector  separations 
are  30  pm  in  the  horizontal  direction  and  18  pm  in  the  vertical 
direction,  so  the  speckle  patterns  in  the  Fourier  plane  were  somewhat 
oversampled.  The  video  digitizer  samples  each  horizontal  line  at  a 
spacing  of  21.3  /im. 
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For  the  speckle  in  the  image  plane,  the  required  sample  spacings 
are  [5.7] : 


Ax 


(5-13) 


where  f0  is  the  focal  length  of  the  lens  L~  and  D  and  D  are  the 
aperture  dimensions.  Since  ^  =  189  mm,  and,  for  the  small  aperture, 
Du  ^  Dv  =  1.25  mm,  the  required  spacing  is  39  p m.  For  the  large 
aperture,  Du  =  -  5  mm,  and  the  required  detector  spacing  is  10  pm . 
Note  that  this  spacing  is  smaller  than  the  pixel  spacing  of  the 
detector  array  used.  Because  this  high-resolution  image  is  only  used 
for  comparison  purposes,  sampling  it  at  less  than  the  Nyquist  rate  has 
no  effect  on  the  success  of  the  iterative  algorithm. 


An  optical  Fourier  transformation  exists  between  the  optical  fields 
in  planes  P^  and  P 2  of  Figure  5-1  while  a  digital  Fourier 
transformation  (e.g.,  an  FFT)  will  be  used  when  processing  the  data 
collected  in  these  two  planes.  The  scaling  between  these  two  data  sets 
can  be  determined  as  follows  [5.1],  The  Fourier  intensity  data  is 
collected  as  N  by  M  values,  electronically  sampled  by  the  video 
digitizer  at  spacings  Au  =  21.3  pm  and  Av  =  18  pm.  When  this  data  is 
FFTed  to  the  image  domain,  the  resulting  data  will  be  again  N  by  M 
samples  with  spacings,  in  the  image  plane  P of  AxQ  =  Xf^/NAu  and  AyQ 
=  Xf^/MAv.  The  actual  low  and  high  resolution  image  data  is,  again, 
electronically  sampled  by  the  video  digitizer  at  spacings  Ax  =  21.3 
and  Ay  =  18  pm,  so  the  data  must  be  rescaled  to  be  sampled  at  Axq  and 
AyQ.  The  scaling  factors  are 

Ax  Xf„ 

AF  =  TZTE  =  °*84 
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and 


Av  Xf~ 

J  o  _  2 

Ay  M  Ay  Av 


(5-14) 


for  N  =  M  =  256. 

For  phase  retrieval,  the  image  to  be  reconstructed  must  occupy  no 
more  than  N/2  by  M/2  samples.  The  magnification  of  the  image  with 
respect  to  the  object  is  f 2/f i  for  the  setup  shown  in  Fig.  5-1.  For  an 
object  of  extent  dx  in  the  x-direction,  the  image  extent  is  therefore 
f9d  If..  Using  Eq.  (5-12)  with  equality  for  Au,  we  find  that  (N/2)Ax 
=  Xf2/2Au  =  f2cJx/f1  verifying  that  proper  sampling  has  been  achieved. 

5.6  EXPERIMENTAL  RESULTS 

In  one  series  of  experiments,  images  were  reconstructed  using  only 
Fourier  intensity  data.  The  result  of  one  of  these  experiments  is 
shown  in  Figure  5-2.  In  part  (a),  the  square  root  of  the  detected 
Fourier  intensity  data  is  shown.  This  256  by  256  array  of  8  bit  data 
is  the  result  of  averaging  256  frames  of  digitized  CCD  data.  From  this 
data,  the  autocorrelation  of  the  desired  complex-valued  image  was 
computed  via  a  Fourier  transformation.  In  part  (b) ,  the  estimate  of 
the  upper  bound  on  the  image  support  obtained  by  applying  the  triple 
intersection  method  to  the  thresholded  autocorrelation  is  given  [5.5]. 
Part  (c)  shows  the  intensity  of  the  image  reconstructed  from  the 
Fourier  intensity  data  of  part  (a)  and  the  support  of  part  (b) .  For 
comparison,  an  optically  formed  image  is  shown  in  part  (d).  It  can  be 
seen  that  additional  details  such  as  the  dark  stripe  on  the  rectangle 
and  the  shape  of  the  triangle  have  been  reconstructed  correctly 
although  at  lower  than  actual  contrast.  Because  this  image  data  was 
collected  at  a  different  time  than  the  Fourier  intensity  data,  with  a 
slightly  different  aperture,  nd  had  to  be  resampled  for  display,  the 
fine  details,  such  as  the  speckles,  are  not  expected  to  agree  with 
those  in  part  (c).  In  a  second  experiment,  images  were  not 
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Figure  5-2.  Image  reconstruction  from  experimental  data,  (a)  Fourier 
modulus,  (b)  image  support  from  triple  intersection 
method,  (c)  reconstructed  image  intensity,  (d)  approximate 
conventional  image  intensity  for  comparison  to  (c) . 


resolution  image  data  using  the  method  described  in  Chapter  4.  Further 


work  to  improve  data  collection  procedures  for  this  case  is  required. 
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6 

SYSTEM  ANALYSIS 


6.1  INTRODUCTION 

A  larg=  segment  of  this  research  program  was  devoted  to  the 
specification  of  the  hardware  and  operational  characteri sti cs  required 
by  a  deployed  system  for  image  reconstruction  and  target  parameter 
estimation.  In  this  section  of  the  report,  the  results  of  this  system 
analysis  study  are  summarized.  Section  6.2  contains  a  summary  of  the 
results  pertaining  to  parameter  estimation  and  Section  6.3  summarizes 
the  results  for  image  reconstruction. 

6.2  SYSTEM  ANALYSIS  FOR  PARAMETER  ESTIMATION 

In  this  section,  several  factors  that  establish  limits  on  the 
ability  to  measure  target  rotation  rates  using  speckle  are  considered. 
After  an  introduction,  source  requirements,  object  material  properties 
and  receiver  requirements  are  discussed  separately. 

6.2.1  Principles  of  Rotation  Measurement 

Let  us  consider  some  of  the  principles  upon  which  rotation 
measurement  using  speckle  is  based  and  make  a  comparison  to 
conventional  coherent  Doppler  radar.  The  premise  for  this  work  is  that 
various  components  of  an  SDI  target  cluster  rotate  at  different  rates 
and  that  measurement  of  the  rotation  rates  within  a  target  cluster  will 
allow  one  to  identify  the  objects  present.  The  use  of  speckle  to 
measure  rotation  rate  can  be  regarded  as  a  Doppler  technique.  Rather 
than  interfering  coherent  radiation  reflected  from  the  object  with  a 
local  oscillator,  as  in  heterodyning  and  homodyning,  the  speckle 
techniques  are  autodyne  or  self-Doppler  methods.  These  methods  are 
based  on  the  interference  of  light  scattered  from  different  regions  of 
the  object  and  uo  not  employ  a  local  oscillator.  Researchers  have 
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shown  that  even  though  autodyning  avoids  the  comp !  .  i '  ,  isscciated  wit": 
-aving  a  local  oscillator,  it  has  the  same  sens : !  ;  *  v  Heterodyning 

'or  measuring  Doppler  spread  [6.1]. 

Observation  of  the  speckle  pattern  produced  by  illuminating  a 
rough,  rotating  object  with  a  laser  beam  reveals  not  only  that  the 
intensity  at  a  point  in  the  speckle  pattern  dees  exhibit  Doppler 
seating  when  the  object  rotates,  but  also  that  the  beating  is  impart ec 
by  :  bodily  translation  of  the  speckle  pattern.  Moreover,  trie  speck,  e 
translation  is  proportional  to  the  object  rotation.  This  reservation 
suggests  that  rather  than  measuring  rotation  simply  by  temporal 
processing  of  the  intensity  detected  at  a  point  as  with  standard 
autodyning,  one  can  combine  spatial  with  temporal  processing  to  extract 
the  rot  at i on  information.  For  example,  if  the  speckle  pattern  is 
recorded  twice  by  a  2-D  detector  array  with  the  time  interval  between 
recordings  equal  to  At,  the  recordings  can  be  spatially  correlated  to 
extract  the  speckle  translation,  Ax,  which  in  turn  gives  the  object 
angular  rotation  rate,  fl  (in  radi ans/sec) ,  via 


q 


J_  Ax 
2R  At 


where  R  is  the  range. 


(6-1) 


\ 

( 


One  other  difference  between  autodyning  and  heterodyning  worth 
noting  is  that  heterodyning  is  intended  for  measurement  of  relative 
motion  between  the  source  and  target  while  autodyning  is  not  sensitive 
to  their  relative  motion.  Hence,  heterodyning  is  most  sensitive  to 
longitudinal  target  motion  while  autodyning  concentrates  on  differences 
of  motion  within  a  target  or  target  cluster.  For  this  reason  the 
frequency  content  of  the  heterodyne  signal  from  an  object  is  biased 
significantly  by  relative  motion  which  can  complicate  the  extraction  of 
rotation  information. 
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Another  quantity  important  to  the  discussion  of  system  capabilities 
is  the  time  that  it  takes  a  speckle  to  translate  over  a  point  in  the 
detection  plane.  This  time  is  the  inverse  of  the  required  detector 
bandwidth.  With  the  speckle  size  given  by  XR/d,  where  X  is  the 
illumination  wavelength,  R  is  the  object  range  and  d  is  the  object 
diameter,  it  follows  that  the  bandwidth  required  of  each  detector  to 
ensure  Nyquist  sampling  is  given  by 

B  =  4fld/X  (6-2) 

The  system  used  to  conduct  these  rotation  measurements  could 
actually  be  constructed  to  operate  in  a  number  of  possible  modes.  We 
believe  that  the  optimal  system  would  have  a  laser  source  that  operates 
in  a  long  pulse  mode  and  a  2-D  detector  array  for  which  each  element 
has  bandwidth  given  by  Eq.  (6-2).  As  the  object  rotates  the  speckle 
pattern  is  recorded  as  a  function  of  time.  To  compute  the  speckle 
translation,  speckle  data  collected  at  a  minimum  of  two  different  times 
is  required.  The  temporal  separation,  At,  of  this  data  should  be  set 
so  that  the  speckle  translation  is  detectable. 

6.2.2  Source  Requirements 

To  accomplish  target  rotation  measurement  for  the  SDI  mission,  the 
laser  source  must  have  sufficient  output  power  and  coherence  length. 
The  coherence  length  must  be  twice  the  depth  of  the  object  so  that,  at 
a  given  time,  the  light  scattered  from  the  entire  target  adds 
coherently  at  the  detector.  If  this  condition  is  met,  the  speckle 
contrast  is  unity;  otherwise,  speckle  with  reduced  contrast  results  and 
measurement  of  speckle  translation  degrades. 

The  laser  pulse  must  be  long  enough  so  that  the  speckle  translates 
a  detectable  amount.  The  pulse  length  must  also  be  long  enough  so  that 
the  light  at  the  detector  originates  from  the  entire  depth,  AR,  of  the 
object.  If  the  pulse  is  too  short,  only  a  slice  of  the  object  is 


illuminated  at  a  time.  With  each  slice  giving  an  independent  speckle 
pattern,  the  resultant  pattern  integrated  over  the  entire  pulse  has 
reduced  contrast  which  decreases  the  ability  to  measure  rotation.  Even 
when  the  pulse  is  sufficiently  long  to  illuminate  the  entire  depth  at 
once,  the  leading  part  of  the  pulse  (of  length  AR)  and  the  trailing 
part  (of  length  AR)  are  not  useful  because  the  entire  object  is  not 
illuminated.  These  regions  of  the  pulse  act  to  reduce  the  speckle 
contrast  and  reduce  measurement  accuracy.  One  can  define  a  pulse 
utilization  parameter  to  quantify  the  fraction  of  the  pulse  that  is 
useful 

Pulse  Utilization  =  1  -  2AR/cr  (6-3) 

where  c  is  the  speed  of  light  and  r  is  the  pulse  length.  The  value 
given  by  Eq.  (6-3)  should  be  as  close  to  unity  as  attainable,  that  is, 
long  pulses  with  r  »  2AR/c  should  be  used. 

If  the  value  of  Eq.  (6-3)  is  positive,  which  signifies  that  for 
some  fraction  of  time  the  entire  object  is  illuminated,  but  the  value 
is  not  close  to  unity  because  the  pulse  is  short,  one  could  employ 
range  gating  at  the  receiver  to  reject  the  leading  and  trailing  edges 
of  the  pulse.  Such  range  gating,  however,  implies  having  detection 
electronics  that  are  more  complex,  and  of  higher  bandwidth,  than  those 
required  for  long  pulse  operation. 

In  summary,  it  is  expected  that  system  complexity  is  minimized  by 
operation  with  long  pulses  for  which  r  »  2AR/c.  Furthermore,  the 
pulse  length  should  be  long  enough  so  that  the  desired  speckle 
translation  is  observed  within  one  pulse.  Except  for  very  slow 
rotation  rates,  this  condition  will  normally  be  met  by  r  »  2AR/c. 
(For  slow  rotation  rates,  two  pulses  can  be  used  to  collect  data  at  the 
required  two  different  times.)  For  this  system  the  bandwidth  for  each 
detector  in  the  receiver  is  given  by  Eq.  (6-2). 


6.2.3  Object  Material  Properties 

To  perform  rotation  measurements  using  speckle  it  is  assumed  that 
the  object  is  diffusely  reflective  so  that  the  speckle  pattern  is  fully 
developed  and  has  unit  contrast.  The  degree  of  diffuseness  is 
characterized  by  the  scattering  cross-section  as  a  function  of  the 
angle  between  the  object's  local  surface  normal  and  the  illumination 
angle.  As  the  degree  of  diffuseness  decreases,  or  the  angular  extent 
of  the  scattering  cross-section  becomes  narrower,  the  principle  effect 
is  that  the  speckle  size  increases  due  to  the  smaller  effective  object 
size.  The  increase  in  speckle  size  reduces  the  SNR  for  rotation 
measurement  because  the  SNR  for  rotation  measurement  is  proportional  to 
the  square  root  of  the  total  number  of  speckles  in  the  detection 
aperture.  Although  the  SNR  is  reduced,  the  bandwidth  requirement  given 
by  Eq.  (6-2)  is  relaxed  because  of  the  increased  speckle  size. 

6.2.4  Receiver  Requirements 

The  receiver  described  in  Section  6.2.1  for  making  rotation 
measurements  using  speckle  is  a  2-D  detector  array.  The  size  of  the 
array  should  be  large  so  that  it  encompasses  a  large  number  of 
independent  speckles  since  the  SNR  is  given  by  Eq.  (2-17): 


SNR 


_ 4MK _ 

20  +  12/ (N<n»  +  2/(N<n>2) 


1/2 


(6-4) 


where  M  is  the  number  of  independent  speckles  encompassed  by  the 
detector  array,  v  is  the  number  of  independent  frame  pairs  used  to 
measure  the  correlation  from  which  rotation  rate  is  determined,  N  is 
the  number  of  detection  channels  per  speckle  and  <n>  is  the  mean  number 
of  detected  photons  per  detection  channel.  The  number  of  detected 
photons  per  speckle  is  equal  to  N<n>.  The  detector  spacing  should  be 
set  according  to  the  finest  expected  Nyquist  sampling  frequency  to 
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avoid  aliasing,  that  is,  the  spacing  shoo'd  be  less  than  or  equal  to 
\R.2d.  It  must  be  noted,  however,  that  o.-ersamp'  iny  results  in  a  lower 
SNR  when  operating  at  low-light  levels  s :  r  p>  ; .  (2-13): 

SNR  =  [2M  K/N] 1  Z  \  ■  (6-5) 

and  the  optimum  value  of  N  is  the  mi  - 1  value,  two,  allowed  by 
Nyquist  sampling,  when  N<n>  is  held  constant. 

5.3  SYSTEM  ANALYSIS  FOR  IMAGE  RECONSTRUCTION 


In  this  section  the  results  of  tne  s/stem  study  for  image 
reconstruction  are  given.  Section  6.3.1  devoted  to  discussion  of 
the  source  and  sensor  requirements  and  Section  6.3.2  contains  a 
discussion  of  simulation  experiments  for  -  na  ra-  t  zat  i  on  of  image 
reconstruct  ion  error  in  the  presence  of  noise.  Section  6.3.3  contains 
a  brief  description  of  the  computational  requirements  for  Image 
reconstruction. 

6.3.1  Source  and  Sensor  Requirements 

The  basic  hardware  elements  for  performing  image  reconstruction  are 
a  pulsed  laser  source,  a  2-0  array  of  light-bucket  type  detectors  and  a 
low-resolution  imager.  As  detailed  in  Section  4,  the  detector  array 
provides  Object  Fourier  magnitude,  or  pupil-plane  information,  while 
the  low  resolution  imager  picks  off  a  portion  of  the  Fourier  magnitude 
data  using  a  beamsplitter  and  provides  a  low-resolution  diffraction- 
limited  image.  This  low-resolution  image  serves  two  purposes;  first, 
application  of  the  Gerchberg-Saxton  phase  retrieval  algori thm  yields 
the  phase  of  the  optical  field  over  the  region  of  the  pupil  plane 
corresponding  to  the  low- resolution  image.  S°cond,  the  low-resolution 
image  serves  as  a  support  constraint  for  reconstruction  of  the  high- 
resolution  image  using  the  expanding  Fourier  modulus  method  described 
in  Section  4. 
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The  source  and  detector  requirements  are  basically  the  same  as 
those  for  parameter  estimation  discussed  in  Section  6.2,  except  that 
for  good  image  reconstruction  performance  (yet  to  be  quantified)  it  is 
required  that  one  have  a  higher  power  laser  (to  increase  the  SNR  of  the 
data)  than  is  required  for  parameter  estimation.  For  slowly  rotating 
objects  or  non-rotating  objects  it  may  be  possible  to  integrate  returns 
from  several  pulses. 

6.3.2  Simulation  Experiments 

In  this  section  the  results  of  simulations  conducted  to 
cnaracterize  the  performance  of  the  image  reconstruction  algorithm  in 
the  presence  of  noise  are  presented.  The  first  simulation  is  a 
comparison  of  image  quality  obtained  using  phase  retrieval  (under  ideal 
circumstances)  to  conventional  diffraction-limited  imaging.  The  object 
used  was  a  simulated  PBV-RV  complex-valued  object.  For  phase 
retrieval,  the  pupil-plane  Fourier  magnitude  was  computed  with  Poisson 
statistics  to  simulate  low  light  level  noise.  An  image  was  then 
reconstructed  using  the  iterative  transform  algorithm  with  the  exact 
image  domain  support  constraint,  and  the  normalized  rms  error  in  the 
reconstructed  image  intensity  was  computed.  For  the  conventional 
diffraction-limited  image,  Poisson  statistics  were  applied  to  the  image 
and  the  resulting  image-domain  error  was  computed. 

Figure  6-1  compares  image  quality  for  phase  retrieval  and 
conventional  imaging  as  a  function  of  the  number  of  photons  per 
speckle.  For  high  light  levels  (>6000  photons  per  speckle)  both 
methods  produce  very  high  quality  images.  At  the  light  level  of 
roughly  1000  photons  per  speckle  the  error  in  the  phase  retrieval 
image  begins  to  increase  markedly.  It  is  important  to  note  that,  even 
at  the  light  level  of  600  photons  per  speckle,  the  error  of  0.30 
associated  with  phase  retrieval  corresponds  to  an  image  that  does 
contain  recognizable  features.  In  practice,  the  support  will  not  be 
known  exactly,  so  imaging  performance  could  be  poorer  than  for  these 
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experiments.  Thus  Figure  6-1  should  not  be  used  to  establish  the  lower 
limit  on  light  level  when  using  phase  retrieval.  In  future  work  a 
thorough  characteri zation  of  the  algorithm  at  low  light  levels  will  be 
conducted  to  establish  the  lower  limit.  This  study  will  include 
analysis  of  the  sensor  configuration,  the  estimation  of  the  support 
constraint,  and  the  specific  form  of  the  phase  retrieval  algorithm  to 
obtain  an  accurate  estimate  for  the  lowest  allowable  light  level. 

Another  series  of  simulations  was  conducted  to  determine 
reconstructed  image  quality  as  a  function  of  noise  in  the  Fourier 
magnitude  data.  The  exact  image  domain  support  constraint  was  used  and 
Gaussian  noise  was  added  to  the  Fourier  intensity.  The  results  of  this 
simulation  are  shown  in  Figure  6-2. 

6.3.3  Computati onal  Requirements 

The  dominant  computational  burden  in  iterative  phase  retrieval 

algorithms  is  the  2-D  FFT  operation;  two  2-D  FFT's  are  required  per 

iteration.  To  reconstruct  an  image  it  thus  follows  that  the  number  of 

real  additions  and  multiplications  required  is  of  the  order 

20MN  / 1 og^N ,  where  M  is  the  number  of  iterations  required  and  N  is  the 

number  of  pixels  in  the  reconstructed  image.  For  N  =  64  and  M  =  100  it 

0 

follows  that  50  x  10  real  additions  and  multiplications  are  required. 

To  assess  the  ability  of  advanced  hardware  to  accommodate  phase 
retrieval  we  considered  VHSIC  technology  and  in  particular,  the 
Westinghouse  Pipelined  Arithmetic  Unit  (PLAU).  The  PLAU  is  capable  of 
performing  40  x  10  operations  per  second.  For  the  above  example  of  N 
=  64  and  M  =  100  it  follows  that  a  single  PLAU  is  capable  of 
reconstructing  an  image  in  about  1.25  sec.  Future  advances  in 
processor  technology  including  highly  parallel  computing  architectures 
can  be  expected  to  lower  the  time  required  for  reconstruction 
consi derably . 
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ERROR  IN  RECONSTRUCTED  IMAGE 


ERROR  IN  FOURIER  MAGNITUDE 


Figure  6-2.  Normalized  rms  error  in  image  reconstructed  by  phase 
retrieval  versus  normalized  error  in  Fourier  magnitude 
data  (due  to  additive  Gaussian  noise  in  Fourier 
i ntensi ty) . 
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APPENDIX 

ERROR  VARIANCE  RELATIONSHIPS 

We  relate  the  variance  of  the  error  of  an  optical  field  (or  the 
Fourier  transform  of  a  complex-valued  image)  to  the  variance  of  its 
phase  error  for  a  zero-mean  Gaussian-distributed  phase  error. 

Let 

G (u)  =  F(u)  exp [i ^e(u) ]  (A-l) 

be  the  aberrated  optical  field,  where  F  is  the  ideal  optical  field  and 

<pe  is  the  phase  error.  Suppose  <j>e  has  point  statistics  that  are 

Gaussian  zero  mean  with  standard  deviation  a..  First  consider  the  case 

9 

without  normalizing  G.  Then  the  variance  of  the  error  (i.e.,  the  mean 
square  error)  of  G(u)  is 

E2  =  A"1  f  I G (u)  -  F(u) 1 2  d2u 

=  A'1  j  I F(u) 1 2  II  -  exp [i ( u) ] 1 2  d2u 
=  A'1  f  I F(u) 1 2  4  sin2(>e(u)/2]  d2u  (A-2) 

where  A  is  the  area  of  integration.  Assuming  the  phase  errors  are 
independent  of  I F ( u ) I ,  and  approximating  the  integral  by  an  ensemble 
average  yields 


E2  =  4  <1 F(u) I  2  sin2  [0e(u)/2]> 

=  4  <IF(u)I2>  <sin2[^(u)/2]>  .  (A-3) 
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Using  the  identity 


I  exp(-px2)  sin2(ax)  dx  =  ^  W  [1  -  exp(-a2/p)]  ,  [p  >  0]  (A-4) 


the  average  over  the  distribution  of  phases  is  given  by 


<sin2[>e(u)/2]>  =  j  sin2(0e/2)  — — -  exp[-02/2a2]  d^ 


=  (1/2)  [l  -  exp[-a2/2)]  . 


(A-5) 


Inserting  this  into  Eg.  (A-3)  yields 


- ^ -  =  2[l  -  exp(-a2/2]l 

<  I  F  ( u )  1 2>  1  *  JJ 


(A-6) 


Note  that  e  -►  2  for  a ,  -»  *  and 

9 


e2  ^  a2  for  aj;  «  1  . 


Next  consider  the  case  of  a  normalized  G,  as  in  Eg.  (4-3): 


(A-7) 


where 


E2  =  A'1  f  la  G(u)  -  F (u) 1 2  d2u 
f  G*(u * )  F(u')  d2u ' 


j  IG(u") 1 


2  d2u" 


(A-8) 


( A- 9 ) 
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<*>  -  - 
<ier> 


<1  Fl^>  <exp(-i^)) 
<IF!2> 


=  <exp(-i^e)>  =  exp(-a^/2)  . 


(A— 10) 


(A— 1 1 ) 


A"1  J  i F ( u) 1 2  lexp(-a2/2)  exp[i^e(u)]  -  II2  d2u 

A"1  J  IF(u)l2  (exp(-a^)  +1-2  exp(-a2/2)  cos  [*e(u)]}  d2u 

< I F (u) I 2>  |exp(-a2)  +1-2  exp(-<72/2)  <cos[^e(u)]>} 

< I F ( u) I 2>  [exp(-a2)  +1-2  exp(-ff2/2)  exp(-a2/2)] 

< I  F ( u )  1 2>  [l  -  exp (-<x2)  ,  (A-12) 


<1  F(u) I 2> 


2  1  -  exp(-<72)  . 


(A- 13) 


Note  that  for  the  normalized  case,  unlike  the  unnormalized  case,  e2  +  1 

2  2  2 

for  a .  +  o#  and,  like  the  unnormalized  case,  e  2  a,  for  a  ,  «  1. 
r  9  <P 

By  Parseval's  theorem  it  can  be  shown  that  the  variance  of  the 
error  in  the  image  domain  is  equal  to  the  variance  of  the  error  in  the 
Fourier  domain. 


Just  as  image  shifts  can  be  taken  out  before  computing  errors  to 
allow  for  the  fact  that  image  shifts  are  unimportant  to  image  quality, 

linear  components  of  the  phase  error  (u)  can  be  taken  out  before 
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computing  ^  or  e  . 


